Profesor: Julio Antonio Soto
Scikit-Learn: machine learning hecho sencillo con Python¶
Ahora que ya tenemos una base del lenguaje Python y conocemos las posibilidades de Numpy y pandas (así como de visualización), ha llegado la hora de aplicar todo lo aprendido al mundo del aprendizaje automático.
Scikit-Learn es, probablemente, la biblioteca de machine learning más extensa que hay. Rápida, estable y en continua evolución, Scikit-Learn (se abrevia como sklearn) es la respuesta evidente a hacer machine learning con Python.
Sklearn está escrita en una mezcla de Python, Cython, C y C++. Utiliza extensamente Numpy y Scipy por debajo, y su integración con nuestras Series y DataFrames de pandas es tan sencilla como te imaginas.
1. Introducción¶
Sklearn es una biblioteca que contiene una gran variedad de modelos escritos, y listos para ser entrenados y utilizados:
- Modelos de clasificación
- Modelos de regresión
- Modelos de clustering
- Modelos de reducción de dimensionalidad
- Utilidades para validación y selección de modelos
- Utilidades para preprocesado de datos y feature engineering
- Una arquitectura de pipelines para hacer workflows complejos de todos los puntos anteriores
A todo el mundo que conozco le encanta el siguiente diagramita que aparece en la página web de sklearn, sobre "cómo elegir tu modelo ideal":
Evidentemente se trata de una ultra-generalización, pero no por ello deja de ser bonita. Y no se puede decir que sea errónea tampoco...
Iremos haciendo referencias a la documentación de sklearn. Solo diré que es probablemente la mejor documentación que he visto, y solo con ella puedes aprender muchísimo sobre machine learning. No, no es broma.
Suficiente. Pongámonos manos a la obra.
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use("ggplot")
2. Cosas importantes antes de ponerse a machacar datos¶
Estarás deseando empezar a lanzar modelos... Pero antes vamos a ver algunas de las opciones que nos ayudarán a hacerlo:
2.A. Creación de datos sintéticos¶
A veces... no tenemos datos. Es lo que hay. Pero eso no tiene que pararnos los pies, sobre todo cuando se trata de aprender. De hecho, la mayoría de I+D en machine learning e inteligencia artificial se hace con datos sintéticos (es decir, generados de forma aleatoria por nosotros); o bien con datasets muy conocidos.
Afortunadamente, scikit-learn trae de serie una serie de funciones que nos permiten generar datos sintéticos de forma muy sencilla.
Por ejemplo: la función make_blobs genera muestras gausianas (es decir, distribuciones normales en dos o más dimensiones) isotrópicamente distribuidas. Vamos a crear tres blobs:
# La biblioteca scikit-learn se llama sklearn
# (para el tema de imports y tal):
from sklearn.datasets import make_blobs
mi_primer_dataset = make_blobs(n_samples=10, # Numero de puntos que creamos
n_features=2, # Numero de atributos (variables)
centers=3, # Numero de centros (numero de gausianos) -> 3 blobs
cluster_std=1.5, # Desviacion tipica de los datos
random_state=2) # Si queremos poner un "seed" al generador de números aleatorios
mi_primer_dataset
Nos devuelve una tupla con dos arrays de Numpy: la primera son las coordenadas de cada punto, y la segunda a qué centro pertenece cada punto; de forma que el punto [ 0.96455381, 0.46894968] pertenece al centro 1, el punto [ -1.21779287, -11.15836353] pertenece al centro 0, y así.
Dibujemos con matplotlib. Preparamos un DataFrame sencillito con los datos (para verlo más claro; aunque podríamos dibujar a pelo con matplotlib y las arrays de Numpy):
# Primero vamos a meter en un DataFrame las coordenadas de cada punto:
dataset_blobs = pd.DataFrame(mi_primer_dataset[0]) # Primer elemento de la tupla
# Y vamos a renombrar las columnas:
dataset_blobs.columns = ["x", "y"]
# Y podemos añadir el centro al que pertenece cada uno al DataFrame.
# ¿Por qué no? Es común llamarlo "label":
dataset_blobs["label"] = mi_primer_dataset[1]
dataset_blobs
# Los colores por defecto de Matplotlib cuando
# hay distintos puntos son feos...
import matplotlib.cm as color_maps
# Ahora podemos usar estos:
# http://matplotlib.org/examples/color/colormaps_reference.html#color-example-code-colormaps-reference-py
dataset_blobs.plot(kind="scatter",
x="x",
y="y",
c="label",
cmap=color_maps.brg,
figsize=(8,6))
pass
Sosísimo. Muy pocos puntos. Vamos a hacer lo mismo, pero aumentando el número de puntos:
mi_primer_dataset = make_blobs(n_samples=250, # 250 puntos será más vistoso.
n_features=2,
centers=3,
cluster_std=0.5,
random_state=2)
dataset_blobs = pd.DataFrame(mi_primer_dataset[0])
dataset_blobs.columns = ["x", "y"]
dataset_blobs["label"] = mi_primer_dataset[1]
dataset_blobs.plot(kind="scatter",
x="x",
y="y",
c="label",
cmap=color_maps.brg,
figsize=(8,6))
pass
Mejor. Hacer un clustering con este dataset sería coser y cantar.
El resto de generadores de datos sintéticos funcionan prácticamente igual que make_blobs. Aquí está la lista completa.
2.B. Datasets precargados¶
Pero siempre podemos preferir tirar de un dataset real. Si no tenemos el nuestro propio, sklearn viene con unos cuantos precargados, así como con utilidades para descargar automágicamente de repositorios online como puede ser mldata.org.
Vamos a cargar un famoso dataset de precios de casas en Boston:
from sklearn.datasets import load_boston
datos_boston = load_boston()
Con load_boston() hemos cargado un diccionario de Python, que tiene las siguientes claves:
datos_boston.keys()
En DESCR supongo que estará la descripción del dataset...
print(datos_boston["DESCR"])
Pues sí. feature_names nos da los nombres de las columnas, data son las variables independientes, y target es el valor dependiente (el que queremos predecir, es decir: el valor de cada casa). Cargarlo todo en un DataFrame es coser y cantar:
# Cargamos las columnas en un DataFrame:
dataset_boston = pd.DataFrame(datos_boston["data"])
# Y ponemos los nombres a las columnas
# (damos por hecho que el orden de las columnas
# es el correcto):
dataset_boston.columns = datos_boston["feature_names"]
# Y añadimos la columna de precios
# (que si leemos más arriba en la descripción,
# podemos ver que quieren que la columna se llame
# MEDV. Pero somos unos rebeldes):
dataset_boston["Precio"] = datos_boston["target"]
dataset_boston[:10]
Distintos datasets precargados pueden venir en distintos formatos; pero siempre utilizando una estructura de dato de Python (probablemente algún tipo de colección) que ya conocemos, así que no deberíamos tener problemas en cargar cualquiera de ellos.
2.C. Separación de datos de training y testing¶
Es común no utilizar todos los datos de los que disponemos para entrenar nuestros modelos. La causa es bien sencilla: los algoritmos de aprendizaje automático tienden a "aprenderse de memoria los datos" con los que entrenan. Cuando entrenamos un modelo, lo más frecuente es que dicho modelo intente minimizar el error cometido dentro de ese dataset de entrenamiento ($E_{in}$, de "dentro" de los datos mostrados al modelo). La gracia del aprendizaje automático es que nuestros modelos sean capaces de generalizar, es decir, que sean capaces de predecir bien en datos a los que no han sido expuestos anteriormente (el error cometido con datos "nuevos" con los que un modelo no ha entrenado se puede escribir como $E_{out}$).
En un mundo ideal los modelos generalizan perfectamente, de forma que $E_{in} \approx E_{out}$, es decir: el error cometido con los datos de entrenamiento es (a ser posible) pequeño, y es igual de pequeño con datos nuevos.
Verificar la diferencia entre $E_{in}$ y $E_{out}$ es muy importante para ver la precisión final que tendrá nuestro modelo. Por eso es buena costumbre reservar parte de los datos y apartarlos del entrenamiento. De forma que:
- Utilizaremos una parte de los datos llamada training set para entrenar los modelos.
- Utilizaremos el resto de los datos (test set) para ver cómo de bueno es nuestro modelo con datos con los que no ha entrenado.
Más adelante entraremos en el tema de la validación.
El caso es que sklearn nos permite dividir nuestro dataset en training y test sets de una forma muy sencilla, con sklearn.model_selection.train_test_split.
Ya que tenemos el dataset de los precios de Boston cargado, vamos a separarlo en training set (80% de los datos) y test set (20% de los datos):
from sklearn.model_selection import train_test_split
datos_spliteados = train_test_split(dataset_boston,
train_size=0.8, # 80% training
test_size=0.2 # 20% testing
)
datos_spliteados es ahora una tupla que contiene los dos DataFrames: el de training y el de testing:
boston_training_set = datos_spliteados[0]
boston_test_set = datos_spliteados[1]
print("Training set (filas, columnas): " + str(boston_training_set.shape))
print("Test set (filas, columnas): " + str(boston_test_set.shape))
2.D. Feature standarization¶
Es frecuente encontrarse con datasets que tienen variables que pertenecen a dominios completamente distintos, y con variopintas magnitudes. Por muy inteligentes que sean nuestros modelos, mezclar peras con manzanas no siempre es ideal. Por ejemplo: si queremos hacer un modelo de scoring de riesgo para clientes, nos podemos encontrar con variables como edad (que probablemente tome valores entro 18 y 100 años) y nivel de ingresos (que estará en euros/miles de euros). Distintas magnitudes y, sobre todo, distintos valores absolutos. En un caso como este es posible que los modelos decidan "dar más importancia" a la variable ingresos; que seguro va a tener unos valores máximos y mínimos más extremos, y probablemente más varianza. Seguro que también va a haber outliers más drásticos...
Una buena solución para esto es lo que se llama feature standarization. Feature standarization es una técnica que sirve para que todas las variables hablen el mismo idioma, y tomen valores relativos similares entre ellas. Además, por lo general, ayuda a que los modelos converjan (o terminen siendo precisos) en menos tiempo e iteraciones.
Hay varias formas de hacer feature standarization. No obstante, la más utilizada es la que hace que todas las variables tengan media $0$ y desviación típica $1$. Y resulta que sklearn tiene una función que hace exactamente esto. En la bibliografía e Internet podrás encontrar sitios donde ponen que: feature standarization = feature demeaning + feature scaling (es decir, hacer dos cosas: quitar media y escalar varianza/desv.típica).
sklearn.preprocessing incluye varias funciones para hacer este tipo de transformaciones iniciales a nuestras variables. Para conseguir que las variables tengan media $0$ y desviación típica $1$ utilizamos sklearn.preprocessing.StandardScaler.
Continuemos con los datos de Boston: vamos a coger el training set, y ver la media y desviación típica de cada variable. Con un .describe() en el DataFrame bastará:
boston_test_set.describe()
Podemos ver que tenemos variopintas medias y desviaciones típicas en cada variable. Vamos a estandarizarlas con StandardScaler:
from sklearn.preprocessing import StandardScaler
# Creamos una instancia de StandardScaler:
scaler = StandardScaler()
# Aplicamos el método .fit(), pasando como
# argumento nuestro dataset:
scaler.fit(boston_training_set)
# Ese .fit() computa la media y desviación
# típica de cada variable, dejando preparado
# el escalador. Para aplicarlo sobre nuestros
# datos, llamamos ahora al método .transform()
# (ha de ser después de .fit(), claro):
boston_training_set_scaled = scaler.transform(boston_training_set)
boston_training_set_scaled debería tener ahora nuestras mismas variables, pero con media $0$ y desviación típica $1$. Vamos a verlo:
# Resulta que scaler.transform() devuelve un array/matriz de Numpy,
# no un DataFrame de pandas. Podemos coger las columnas del
# boston_training_set original, y poner los nombres:
boston_training_set_scaled = pd.DataFrame(boston_training_set_scaled)
boston_training_set_scaled.columns = boston_training_set.columns
boston_training_set_scaled.describe()
Las medias son ahora valores muy cercanos a $0$; y la desviaciones típicas valores muy cercanos a $1$. Por tanto, nuestras variables ahora están estandarizadas.
Nota: con el fin de simplificar el código, en este ejemplo hemos estandarizado también la target variable, variable dependiente o variable a predecir (o sea, el precio de las casas). En un escenario real no vas a querer hacer esto; solo vas a querer estandarizar las variables que no quieres predecir.
Sklearn incluye otros varios transformadores y preprocesadores de datos. Puedes leer más sobre ellos aquí y aquí.
Con conocimientos de Numpy, pandas y sklearn.preprocessing no deberíamos tener problemas para preparar nuestros datasets para el machine learning que va luego.
Sin más dilación, empecemos con los modelos...
3. Clasificación¶
Desde mi experiencia, diría que el 80% de los problemas de machine learning/data science que nos encontramos por ahí son problemas de clasificación. Los seres humanos tendemos a "etiquetar" todo dentro de una serie de grupos controlados, con el fin de hacer las cosas más sencillas para nuestro consciente. De hecho, nos pasamos el día clasificando cosas de forma involuntaria. Gracias a eso, podemos ver dos mesas distintas; pero sabemos que ambas son mesas (clasificamos las dos como mesas; diferenciándolas de otras labels: sillas, camas...).
Los problemas de clasificación son sencillos de entender: tenemos un dataset con una serie de variables independientes o features, y una variable que es la clasificación, clase o label de cada observación. Con esos datos, entrenamos en modelo. Cuando en el futuro nos lleguen únicamente las features, seremos capaces de predecir la label de cada observación.
Los problemas de clasificación se resuelven mediante apendizaje supervisado. Esto quiere decir que tenemos en nuestro dataset de entrenamiento las labels de cada observación, y las utilizamos para "enseñar" al modelo.
3.A. Preparemos un dataset sencillo e intuitivo¶
Vamos a crear un dataset simple para recrear un problema de clasificación en dos dimensiones (¡que podamos dibujar!). Vamos a generar datos sintéticos. Podríamos utilizar make_blobs para este propósito. Pero los desarrolladores de sklearn son muy listos, y saben que la clasificación es el problema más común en machine learning. De forma que han creado una función que genera automáticamente datasets sintéticos para problemas de clasificación: sklearn.datasets.make_classification
from sklearn.datasets import make_classification
dataset = make_classification(n_samples=600, # Número de observaciones
n_features=2, # Número de features
n_redundant=0, # Número de features "redundantes" que añaden ruido
n_classes=2, # Número de distintas clases o distintas labels
random_state=749 # Seed del generador de números aleatorios
)
Al igual que con make_blobs, el resultado es una tupla con dos arrays de Numpy: la primera tiene las features, y la segunda las labels:
# Metemos en un DataFrame las features primero:
dataset_clasificacion = pd.DataFrame(dataset[0])
# Y ahora las labels:
dataset_clasificacion["label"] = dataset[1]
dataset_clasificacion[:5]
Ahora inventémonos una historia para nuestros datos. Supongamos que este dataset muestra en cada fila la observación de los niveles de concentración de dos hormonas masculinas (hormona_a y hormona_b) de un hombre. Y supongamos que la label que acompaña a las medidas de hormonas de cada hombre nos dice si dicho hombre se ha quedado calvo o no (donde $0$ es hombre con pelo, y $1$ es hombre calvo). Vamos a intentar, por tanto, predecir si un hombre se quedará calvo o no a partir de sus niveles hormonales.
Disclaimer: DataHack no se hace responsable del posible mal gusto del ejemplo de clasificación elegido.
# Vamos a renombrar las columnas de features:
dataset_clasificacion.columns = ["hormona_a", "hormona_b", "label"]
# Y sumar 4 a los valores de ambas features
# (porque hay concentraciones negativas
# de hormonas si no, y es raro):
dataset_clasificacion["hormona_a"] = dataset_clasificacion["hormona_a"] + 4
dataset_clasificacion["hormona_b"] = dataset_clasificacion["hormona_b"] + 4
dataset_clasificacion[:5]
dataset_clasificacion.describe()
Dibujemos:
# Queda mejor si pintamos ceros y unos por separado:
dataset_clasificacion_ceros = dataset_clasificacion[dataset_clasificacion["label"] == 0]
dataset_clasificacion_unos = dataset_clasificacion[dataset_clasificacion["label"] == 1]
axes = dataset_clasificacion_ceros.plot(kind="scatter",
x="hormona_a",
y="hormona_b",
label="No calvo",
c="blue", # Color de los puntos
edgecolors="black", # Color del contorno de los puntos
linewidths=1, # Ancho del contorno de los puntos
s=30, # Tamaño de los puntos
figsize=(8,6))
dataset_clasificacion_unos.plot(kind="scatter",
x="hormona_a",
y="hormona_b",
label="Calvo",
c="red",
edgecolors="black",
linewidths=1,
s=30,
ax=axes)
pass
Donde cada punto es un señor con sus dos niveles de concentración de hormonas, y si es calvo o no.
3.B. Funcionamiento de los clasificadores¶
Casi siempre vamos a querer hacer lo mismo:
- Cargar un modelo (hacer el
importcorrespondiente) - Entrenarlo con nuestros datos
- Probar su rendimiento
- Predecir nuevos casos (en este caso, si señores nuevos se quedarán o no calvos)
Pues bien, es muy sencillo. Dado un modelo:
- Hacemos
from sklearn.loquesea import modelo modelo.fit(features, labels)modelo.score(features, labels)modelo.predict(features, labels)
El orden en el que queremos hacer los pasos puede cambiar, pero por lo general será algo así.
3.C. Separamos el dataset en training y testing¶
Pues eso. Igual que en el ejemplo de Boston:
dataset_clasificacion_spliteado = train_test_split(dataset_clasificacion,
train_size=0.8, # 80% training
test_size=0.2, # 20% testing
random_state=22 # seed para el generador de num. al.
)
dataset_clasificacion_training = dataset_clasificacion_spliteado[0]
dataset_clasificacion_testing = dataset_clasificacion_spliteado[1]
Y dibujamos:
plt.figure(figsize=(16,7))
# Training
dataset_clasificacion_training_ceros = dataset_clasificacion_training[dataset_clasificacion_training["label"] == 0]
dataset_clasificacion_training_unos = dataset_clasificacion_training[dataset_clasificacion_training["label"] == 1]
axes1=plt.subplot(1,2,1, title="Training")
dataset_clasificacion_training_ceros.plot(kind="scatter",
x="hormona_a",
y="hormona_b",
label="No calvo",
c="blue",
edgecolors="black",
linewidths=0.5,
s=30,
ax=axes1)
dataset_clasificacion_training_unos.plot(kind="scatter",
x="hormona_a",
y="hormona_b",
label="Calvo",
c="red",
edgecolors="black",
linewidths=0.5,
s=30,
ax=axes1)
# Testing
dataset_clasificacion_testing_ceros = dataset_clasificacion_testing[dataset_clasificacion_testing["label"] == 0]
dataset_clasificacion_testing_unos = dataset_clasificacion_testing[dataset_clasificacion_testing["label"] == 1]
axes2 = plt.subplot(1,2,2, title="Testing")
dataset_clasificacion_testing_ceros.plot(kind="scatter",
x="hormona_a",
y="hormona_b",
label="No calvo",
c="blue",
edgecolors="black",
linewidths=0.5,
s=30,
ax=axes2)
dataset_clasificacion_testing_unos.plot(kind="scatter",
x="hormona_a",
y="hormona_b",
label="Calvo",
c="red",
edgecolors="black",
linewidths=0.5,
s=30,
ax=axes2)
pass
3.D. Entrenamiento de clasificadores¶
No podemos saber qué modelo de clasificación va a ser el que mejor funcione para un determinado dataset. Podemos tener nuestras hipótesis o corazonadas; pero fiarse solo en eso contradice el método científico...
... De forma que tendremos que probar distintos clasificadores, y ver cuál es el que mejor funciona. Estos son algunos de los clasificadores que incluye sklearn (hay más):
| Algoritmo | Módulo | Multiclase | Comentario |
|---|---|---|---|
| Regresión Logística | linear_model.LogisticRegresion |
OvR | Emplea la biblioteca liblinear |
| Árboles de decisión | tree.DecisionTreeClassifier |
Inherente | Clasificador de reglas |
| K Nearest Neighbors | neighbors.KNeighborsClassifier |
Inherente | Clasificador en función de la clase de las muestras más cercanas |
| Naïve Bayes | naive_bayes.GaussianNB |
Inherente | Clasificador bayesiano que asume independencia entre atributos |
| Perceptrón | linear_model.perceptron |
OvR | Clasificador linear de margen |
| SVM (sólo lineal) | svm.LinearSVC |
OvR | Emplea la biblioteca liblinear. Sólo kernel lineal |
| SVM | svm.SVC |
OvO | Emplea la biblioteca libsvm. Todo tipo de kernel |
Vamos a ir probando entonces:
3.D.a. Regresión logística¶
A pesar de su nombre la regresión logística no es un algoritmo de regresión, sino de clasificación. Su formulación matemática se basa en la de la regresión lineal, pero luego aplica una función (función logística) para convertir el resultado de dicha regresión en números entre 0 y 1 (por lo general, aunque existen variantes que lo hacen entre -1 y 1). Gracias a este comportamiento podemos interpretar el resultado de predicción como un resultado probabilístico (las probabilidades siempre están entre 0 y 1). Una regla sencilla es que si el resultado > 0.5, se predice como 1, y si el resultado es < 0.5, se predice como 0.
En prácticamente todos los modelos existen una serie de hiperparámetros y argumentos que tenemos que fijar. La formulación matemática de cada modelo será distinta y, por tanto, los resultados también serán distintos, en función de los hiperparámetros que elijamos. La selección de estos hiperparámetros es un mundo complejo, ya que altera la especificación de cada modelo; incluso puede añadir regularizaciones a los mismos. Esta selección de hiperparámetros permite evitar el sobreajuste de los modelos.
Comencemos por entrenar una regresión logística. Evidentemente, con nuestros datos de entrenamiento. El modelo es sklearn.linear_model.LogisticRegression:
from sklearn.linear_model import LogisticRegression
# Creamos una instancia del modelo,
# con todos los hiperparámetros
# y argumentos por defecto:
lr = LogisticRegression()
# Ahora vamos a entrenarla con el conjunto de training:
lr.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]], # nuestras dos features
y=dataset_clasificacion_training["label"]) # labels
El modelo ha sido entrenado, y vemos que devuelve un montón de opciones e hierparámetros elegidos (que puesto que no hemos elegido nada, ha dejado los que vienen por defecto).
Con el modelo entrenado, vamos a probar a hacer un .predict() sobre los datos de testing, para intentar predecir cuáles están calvos/no calvos. Evidentemente tenemos las labels reales de los datos de testing, y luego ya entraremos en el tema de comprobar rendimiento. De momento queremos ver que con nuevos datos que el modelo no ha visto, nuestra regresión logística predice ceros o unos:
lr.predict(dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
Para cada observación de nuestro conjunto de testing, predice 0 o 1. Pero nosotros tenemos los valores reales...
np.array(dataset_clasificacion_testing["label"])
Mirando un poco por encima, vemos que no acierta el 100% en el conjunto de testing. No pasa nada; ahora lo miramos en más profundidad. De hecho, es tan sencillo como utilizar lr.score(features_testing, labels_testing). Básicamente hace lo que hemos hecho nosotros: hace lr.predict() sobre los datos de testing, y compara los resultados de las labels predichas con las reales:
lr.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"])
Hemos conseguido un número mayor que 0.9. Esto quiere decir que hemos acertado si es 0 o 1 en más del 90% de los casos del dataset de testing (nuestro $E_{out}<10\%$). Esta métrica es conocida como accuracy.
Vale: hemos entrenado la regresión logística y hemos comprobado su accuracy con datos con los que no había entrenado (el conjunto de testing). Si recuerdas la regresión lineal, el resultado del modelo era una función lineal:
$$f(x) = \beta_{0} + \beta_{1}{hormona_a} + \beta_{2}{hormona_b}$$Pues en la regresión logística es exactamente igual, solo que con la aplicación de la función logística al final. Vamos a ver los parámetros $\beta$ estimados por el modelo:
print(lr.coef_) # beta1 y beta2
print(lr.intercept_) # beta0, también llamado bias o intercept
print("f(x) = " + str(lr.intercept_[0]) + " + " + str(lr.coef_[0][0]) + " hormona_a + " + str(lr.coef_[0][1]) + " hormona_b")
Esto sería para la regresión lineal normal. Para la regresión logística, tenemos que aplicar la función logística al final (en los libros la veréis también como función sigmoidal o función sigmoide):
$$\theta(f(x)) = \frac{e^{f(x)}}{1+e^{f(x)}}$$Y así es como conseguimos que nuestro valor predicho $\theta(f(x))$ esté entre 0 y 1.
Te preguntarás si se puede dibujar esto de forma sencilla. La gente de sklearn ha preparado algunas funciones (como esta, esta, esta y esta) para que los clasificadores sean lo más fáciles posible de visualizar. Siéntete libre de utilizar (¡y mejorar!) esta función, que cogí de ellos y alteré ligeramente. Dibuja nuestro dataset, así como las "regiones de decisión" de nuestro clasificador:
def plot_decision_regions(clf, X, y, fig=None, title='', xlabel='', ylabel='', figsize=(8,6)):
"""
Dibuja las regiones de decisión de un clasificador de sklearn.
Toma como argumentos:
- clf: la instancia del clasificador entrenado.
- X: las features del dataset.
- y: las labels del dataset.
- fig: la plt.figure() o axes a usar. Si es None, se creará una figure nueva.
- title: el título del gráfico.
- xlabel: el nombre del eje x.
- ylabel: el nombre del eje y.
- figsize: si fig es None, el tamaño de la figura a crear (ancho, alto)
Devuelve:
Nada. Simplemente dibuja el resultado con Matplotlib.
"""
from matplotlib.colors import ListedColormap
cmap_light = ListedColormap(list(reversed(['#FFAAAA', '#AAFFAA', '#AAAAFF'])))
cmap_bold = ListedColormap(list(reversed(['#FF0000', '#00FF00', '#0000FF'])))
# Plot the decision boundary. For that, we will assign a color to each
# point in the mesh [x_min, m_max]x[y_min, y_max].
h = .05 # step size in the mesh
x_min, x_max = X.values[:, 0].min() - .5, X.values[:, 0].max() + .5
y_min, y_max = X.values[:, 1].min() - .5, X.values[:, 1].max() + .5
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
# Put the result into a color plot
Z = Z.reshape(xx.shape)
if fig is None:
figura = plt.figure(figsize=figsize)
plt.pcolormesh(xx, yy, Z, cmap=cmap_light, alpha=0.4)
plt.scatter(X.values[:, 0], X.values[:, 1], c=y.values,
cmap=cmap_bold, edgecolors='k')
plt.title(title)
plt.xlabel(xlabel)
plt.ylabel(ylabel)
#plt.xticks(())
#plt.yticks(())
plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
else:
fig.pcolormesh(xx, yy, Z, cmap=cmap_light, alpha=0.4)
fig.scatter(X.values[:, 0], X.values[:, 1], c=y.values,
cmap=cmap_bold, edgecolors='k')
fig.set_title(title)
fig.set_xlabel(xlabel)
fig.set_ylabel(ylabel)
fig.set_xticks(())
fig.set_yticks(())
fig.set_xlim(xx.min(), xx.max())
fig.set_ylim(yy.min(), yy.max())
No es la función más bonita ni agradable de leer del mundo. Pero vamos a probarla para dibujar nuestra clasificación en los datos de training y testing:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Regresión logística", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(lr,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(lr,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
Así es como nuestra regresión logística ha aprendido a separar los casos de calvo frente a los de no calvo. Aprendió a hacerlo en el training set, y hemos medido su rendimiento en el testing set. Vemos que la regresión logística separa linealmente los datos (con una línea recta). Puede que si te fijes bien veas unos pequeños escaloncitos en el dibujo de la línea recta, pero es tema de la visualización: la línea de separación es recta en la realidad.
Podemos ver también que en el conjunto de testing no acierta a predecir todos los puntos (hay puntos rojos en la zona donde el modelo decide que esos son azules, y viceversa). Aún así, hemos comprobado con lr.score(features_testing, labels_testing) que este error $E_{out}$ es menor al 10%. El modelo ni siquiera es capaz de clasificar bien todos los puntos en el dataset de entrenamiento... ($E_{in}$ no es 0%).
No obstante, la regresión logística es uno de los clasificadores más utilizados. Es sencillo de entender y de implementar, y es bastante rápido entrenando.
Muy rápido, vamos a probar a tocar alguno de los (hiper)parámetros de la regresión logística, para ver si afecta mucho al resultado. Para ello, creamos una nueva:
lr2 = LogisticRegression(penalty="l2", # Penalización L2 o Ridge
solver="lbfgs", # El algoritmo que se usa para minimizar el error, y entrenar el modelo
max_iter=600) # Número máximo de iteraciones para el entrenamiento
Entrenamos con el conjunto de training:
lr2.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
Comprobamos el rendimiento en el conjunto de testing:
lr2.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"])
Y dibujamos:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Regresión logística (regularización L2, L-BFGS, 600 iteraciones)", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(lr2,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(lr2,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
No ha cambiado mucho la cosa (algo ha mejorado, sí). En problemas de 2 dimensiones (2 variables/features, en este caso hormona_a y hormona_b), las posibles regularizaciones que podamos elegir no suelen tener mucho efecto en la regresión logística (especialmente si los datos no tienen valores extremos o mucho ruido).
Por ello, vamos a seguir con el mismo dataset y probar otros clasificadores.
3.D.b. Árboles¶
Los árboles de decisión son muy sencillos de entender, rápidos de entrenar (cuando es de uno en uno; luego cuando veamos ensembles veremos que no son tan rápidos) y son modelos no lineales, es decir: que la línea que separe los datos no tiene que ser necesariamente recta.
Importemos el modelo y entrenémoslo. El modelo es el sklearn.tree.DecisionTreeClassifier:
from sklearn.tree import DecisionTreeClassifier
arbol1 = DecisionTreeClassifier(criterion="entropy", # Criterio de división (entropía o Gini soportado)
max_depth=5, # Máxima profundidad de árbol
min_samples_split=1, # Número mínimo de puntos en un nodo para hacer split
min_samples_leaf=1 # Número mínimo de puntos para considerar un nodo un nodo-hoja
)
arbol1.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
Entrenado. Vamos a ver qué rendimiento nos da en el conjunto de testing:
arbol1.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"])
Y dibujemos:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Arbol de decisión (prof. max. 5)", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(arbol1,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(arbol1,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
Ahora vemos que las zona de decisión es no lineal (esa forma de "rectangulitos" es muy peculiar de los árboles). No obstante, nuestro rendimiento en el conjunto de testing es peor que en la regresión logística, a pesar de que el modelo es más complejo. ¿Por qué?
Porque estamos sobreajustando: el árbol "se está aprendiendo demasiado" los datos de entrenamiento, y cuando tiene que clasificar los del conjunto de testing, comete más fallos que si separara los datos con una sencilla línea recta. Es el principal inconveniente de los árboles: tienden a sobreajustar bastante. Cuando veamos los ensembles veremos distintas formas de evitar que los árboles sobreajusten.
A la gente le suele gustar "ver" los árboles, para ver las reglas de decisión. Sklearn nos permite visualizarlo fácilmente, utilizando sklearn.tree.export_graphviz. La visualización se exporta en formato DOT, el cual podemos pasar a pdf, imagen png y otros. Vamos a ver cómo pasarlo a png y mostrarlo en el notebook.
Nota: antes de instalar las bibliotecas mostradas en la siguiente celda, es posible que necesites instalar graphviz en tu sistema operativo.
from sklearn.tree import export_graphviz
# También necesitamos la biblioteca pydot si usamos Python 2,
# o pydotplus si usamos Python 3. Estas bibliotecas permiten gestionar
# los archivos de formato DOT. Puedes instalar pydotplus con pip:
import pydotplus
# Además, vamos a usar un truco: en vez de exportar el archivo dot
# como archivo al disco duro, vamos a meterlo en una variable
# utilizando la siguiente biblioteca para escribir el output
# del archivo DOT a un string:
from sklearn.externals.six import StringIO
# Creamos una instancia de StringIO, que es un string
# especial donde vamos a "volcar" el archivo DOT:
archivo_dot_en_string = StringIO()
# Exportamos el gráfico:
export_graphviz(arbol1,
out_file=archivo_dot_en_string,
feature_names=["hormona_a", "hormona_b"],
class_names=["no_calvo","calvo"], # Por orden de menor a mayor. 0 es no_calvo, y 1 es calvo
filled=True, # Con colores
rounded=True # Con esquinas redondeadas
)
# Usamos pydotplus para procesar el archivo DOT:
grafo_del_arbol = pydotplus.graph_from_dot_data(archivo_dot_en_string.getvalue())
# Y ahora lo podemos mostrar en el notebook con IPython.display.Image:
from IPython.display import Image
Image(grafo_del_arbol.create_png(), width=900)
Mucho split para solo dos features... Apesta a sobreajuste. Vamos a lanzar uno con menos profundidad:
arbol2 = DecisionTreeClassifier(criterion="entropy",
max_depth=3,
min_samples_split=5,
min_samples_leaf=5
)
arbol2.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
# Score en testing:
arbol2.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"])
Ahora tiene menor $E_{out}$.
figura = plt.figure(figsize=(16,7))
figura.suptitle("Arbol de decisión (prof. max. 3)", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(arbol2,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(arbol2,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
archivo_dot_en_string = StringIO()
export_graphviz(arbol2,
out_file=archivo_dot_en_string,
feature_names=["hormona_a", "hormona_b"],
class_names=["no_calvo","calvo"],
filled=True,
rounded=True
)
grafo_del_arbol = pydotplus.graph_from_dot_data(archivo_dot_en_string.getvalue())
Image(grafo_del_arbol.create_png(), width=900)
Mucho mejor. Pasemos al siguiente clasificador.
3.D.c. Nearest Neighbors¶
Nearest Neighbors (o K-Nearest Neighbors o KNN) es también un modelo muy sencillo de entender. También es no lineal (es decir, que genera regiones de decisión no lineales), y en algunos casos da un gran balance entre sencillez de modelo y resultados. Se basa en un principio muy básico: si conocemos el label de un determinado punto, los puntos cercanos deberían tener el mismo label; y esa será su predicción. Dependiendo de cuántos vecinos contemplemos para un punto a clasificar, el resultado será distinto.
Por ejemplo: supongamos 5-NN, o contemplar los 5 puntos conocidos más cercanos al que queremos clasificar. Si tenemos ese punto a clasificar rodeado de tres calvos y dos no-calvos, se predicirá que el punto a clasificar es calvo (porque de sus 5 vecinos más cercanos, la mayoría (tres) son calvos).
Ahora, para el mismo punto a clasificar, contemplamos 7-NN; y resulta que tres conocidos son calvos (los mismos 3 de antes), pero los otros cuatro conocidos más cercanos son no-calvos (los dos de antes, más otros dos). Ahora, el punto será predicho como no-calvo.
Nearest Neighbors es un modelo no lineal, como veremos en los plots de más abajo.
Sklearn implementa el clasificador KNN a través de sklearn.neighbors.KNeighborsClassifier. Probémoslo:
from sklearn.neighbors import KNeighborsClassifier
knn1 = KNeighborsClassifier(n_neighbors=5, # Número de vecinos contemplados
metric="euclidean", # Distancia utilizada para ver cuáles son los más cercanos
algorithm="brute" # Algoritmo de búsqueda de vecinos (por fuerza bruta)
)
knn1.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
Rendimiento en el conjunto testing:
knn1.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"])
Plots:
figura = plt.figure(figsize=(16,7))
figura.suptitle("5-NN", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(knn1,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(knn1,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
No está mal. Probemos con 1-NN (es decir, contemplar solo el punto más cercano):
knn2 = KNeighborsClassifier(n_neighbors=1, # Número de vecinos contemplados
metric="euclidean", # Distancia utilizada para ver cuáles son los más cercanos
algorithm="brute" # Algoritmo de búsqueda de vecinos (por fuerza bruta)
)
knn2.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
print("Bien clasificados en testing: " + str(knn2.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"])))
# Plot:
figura = plt.figure(figsize=(16,7))
figura.suptitle("1-NN", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(knn2,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(knn2,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
Mayor $E_{out}$, y sobreajuste obvio.
El resultado de K-NN puede cambiar ampliamente en función de la métrica de distancia que utilicemos. La euclídea funciona bien para espacios dimensionales pequeños (este por ejemplo es 2-D); pero cuando tenemos muchas features, es mejor utilizar otra de las métricas de distancia que nos ofrece sklearn.
3.D.d. Naïve Bayes¶
Naïve Bayes son una familia de modelos puramente probabilísticos que son bastante utilizados en clasificación. Como el propio nombre indica, están fundamentados en el teorema de Bayes. Dada una label (en nuestro caso calvo/no-calvo) y una serie de features (en nuestro caso hormona_a y hormona_b), Bayes dice que la relación es:
Dado esto, Naïve Bayes asume que no existe ninguna relación entre la concentración de la hormona a y la de la hormona b. Esto es, que son probabilísticamente independientes (su probabilidad condicionada es la misma que sus probabilidades independientes). Esto es lo que hace Naïve (o ingenuo) al modelo: en la realidad, sí suelen existir relaciones entre las features; sean conocidas o no.
Dentro de sklearn.naive_bayes tenemos varios modelos, en función de cómo asumimos que es la verosimilitud de las features o evidencias, es decir, cómo creemos que es $P(x_i \mid y)$.
Si bien muchos de los otros modelos de machine learning están basados en mayor o menor medida en verosimilitud y modelos probabilísticos, Naïve Bayes está puramente basado en ello. El principal problema de los modelos probabilísticos "puros" es que nos obligan a hacer asunciones sobre cuál es la distribución de probabilidad de las features; y todo lo basado en asunciones, puede ser muy erróneo...
Probemos pues. Vamos a usar un GaussianNB:
from sklearn.naive_bayes import GaussianNB
nb = GaussianNB() # GaussianNB pasa del tema argumentos y parámetros de regularización...
nb.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
print("Bien clasificados en testing: %.15f" % nb.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"]))
# Plot:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Naïve Bayes", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(nb,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(nb,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
Dependiendo de la verosimilitud de las features (en este caso gausiana, de ahí lo de GaussianNB), Naïve Bayes será lineal o no lineal. Por lo general, nos encontraremos que el Naïve Bayes es un clasificador no lineal (como en el ejemplo que acabamos de hacer, aunque no se aprecia demasiado).
3.D.e. Otros clasificadores: SVMs, redes y demás¶
Sklearn incluye otros muchos clasificadores, y sobre todo, quasi-infinitas variaciones de los ya vistos (gracias a la cantidad de hiperparámetros que podemos tocar para cada modelo). Por último, vamos a probar los SVMs (Support Vector Machines), los cuales se verán en profundidad en el módulo de algoritmos avanzados.
Siendo perfectamente discutible: desde mi experiencia, los SVMs son los mejores clasificadores que hay (más que las redes neuronales, ensembles de árboles, y demás), y la matemática detrás de su formulación es una belleza.
Probemos rápidamente un SVM:
from sklearn.svm import SVC
svm = SVC(kernel="rbf") # Kernel de función base radial (ya lo verás)
svm.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
print("Bien clasificados en testing: %.16f" % svm.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"]))
# Plot:
figura = plt.figure(figsize=(16,7))
figura.suptitle("SVM con kernel RBF", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(svm,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(svm,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
Desde la versión 0.18 (28 de septiembre de 2016) Sklearn incluye tres tipos de redes:
La RBM (Restricted Boltzmann machines), que no sirven para clasificación como tal, en su formulación pura (se utilizan sobre todo como autoencoders). Las RBMs son muy buenas aprendiendo distribuciones de probabilidad. Este post de los señores de DeepLearning4J explica muy bien cómo funciona una RBM. La RBM ya estaba disponible en versiones anteriores de sklearn. Las otras dos son nuevas de la versión 0.18.
MLPClassifier implementa la versión clasificador del MLP (Perceptrón Multicapa). El perceptrón multicapa es la red neuronal más sencilla que hay, pero aún así es la más utilizada para datasets convencionales (los que no son de imágenes, texto, secuencias, vídeos... Sino variables "de negocio"). Podemos entender el perceptrón multicapa como una sucesión de capas de regresiones logísticas (más o menos).
MLPRegressor implementa el MLP pero como regresor en vez de clasificador (de forma que la última capa es como las regresiones lineales).
En general, el diseño de redes neuronales requiere tocar muchos más hiperparámetros que el resto de modelos que hemos visto. De ahí que existan bibliotecas específicas para el diseño y entrenamiento de redes neuronales... Y la mayoría de ellas están escritas para Python. Destaco:
- Theano: la de los campeones. Llamarla solo como una biblioteca para hacer redes es casi insultar a Theano. Es una biblioteca capaz de manejar tensores (piensa en arrays de Numpy) de forma simbólica, para delegar la ejecución en tarjetas gráficas o GPUs, que son mucho más rápidas para entrenar redes. Muy manual y complicada: una red sencilla son unas 1.500 líneas de código fácilmente.
- TensorFlow: hecha principalmente por Google, es más sencilla de utilizar. No es mi favorita en absoluto (soy más the Theano), pero va en gustos. Cuando se liberó hace aproximadamente un año era muy pobre, tanto en rendimiento como en funcionalidad. Pero va mejorando, y actualmente existe una biblioteca (también parte del proyecto) llamada TFLearn que es similar a Keras. Independientemente de que la acabes usando o no, existe mucha documentación sobre ella asociada a aprender sobre redes.
- Keras: probablemente la más recomendable para empezar (siempre y cuando ya sepas de redes, y solo te quede por conocer una buena bilioteca para empezar a entrenarlas). Keras es una capa de abstracción por encima tanto de Theano como de TensorFlow. Esto quiere decir: diseñas redes de forma relativamente sencilla con Keras, y se ejecutan utilizando los motores bien the Theano, bien de TensorFlow (recomiendo encarecidamente utilizar el de Theano a día de hoy). Está teniendo un crecimiento brutal en los últimos meses, y es la favorita de muchos.
Ejemplo del MLP de sklearn 0.18 en acción (me tocó usarlo en el curro para una ppt...). Es curioso ver cómo va mejorando el ajuste con las iteraciones:
Vamos a entrenar también un MLP con nuestro dataset (aunque el efecto es menos vistoso, básicamente porque nuestro dataset se puede separar casi con una línea recta):
from sklearn.neural_network import MLPClassifier
# Los hiperparámetros de las redes son difíciles de entender sin haber estudiado bien las redes.
# En el módulo de algoritmos avanzados veréis cómo funcionan estos parámetros:
mlp = MLPClassifier(hidden_layer_sizes=(15,15), # Número de neuronas en las capas ocultas
activation="logistic", # En cada neurona se aplica la función logística
solver="sgd", # Descenso gradiente estocástico
batch_size=10, # Cuán "estocástico" es el descenso gradiente (10 observaciones por "epoch")
alpha=0.0,
learning_rate="adaptive",
learning_rate_init=0.015,
momentum=0.8,
nesterovs_momentum=True,
max_iter=500)
mlp.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
print("Bien clasificados en testing: %.16f" % mlp.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"]))
# Plot:
figura = plt.figure(figsize=(16,7))
figura.suptitle("MLP (dos capas ocultas de 15 neuronas),\nactivación logística/sigmoide", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(mlp,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(mlp,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
Ahora continuaremos con los clasificadores, pero de otra manera. Vamos a empezar con los ensembles.
4. Ensemble Learning¶
El principio del Ensemble Learning es sencillo: la unión hace la fuerza. Consiste en combinar varios modelos en uno solo, con el fin de que el modelo final sea más robusto y generalice mejor.
La intuición es la siguiente: mediante la aplicación de varios modelos a la vez, los fallos de uno de los modelos son compensados por los aciertos de otros de los modelos sobre los mismos puntos. Los ensembles son muy utilizados en la práctica, y la mayoría de los Kaggle se ganan mediante la generación de ensembles de modelos.
Sklearn implementa varias técnicas de ensembles, que se dividen básicamente en dos grupos:
- Bagging, donde se lanzan varios modelos a la vez (normalmente cada uno de ellos con subconjunto del training set). Las predicciones finales se basan en el consenso al que llegan los modelos lanzados: se puede poner a los modelos a votar en el caso de clasificadores, hacer la media de la predicción en el caso de regresores...
- Boosting, donde primero se lanza un modelo. Luego, se lanza otro utilizando parte de las observaciones bien clasificadas por el primero, y las observaciones mal clasificadas por éste (o modificaciones sobre dichos datos). Y así sucesivamente. Se podría decir que el segundo modelo intenta corregir los fallos del primero, el tercero intenta corregir los del segundo, etcétera.
sklearn.ensemble es el lugar donde se encuentran los ensembles. Algunos de ellos son los siguientes:
| Técnica | Módulo | Comentario |
|---|---|---|
| Bagging | ensemble.BaggingClassifier |
Media/votación de múltiples ejecuciones del estimador sobre diferentes subconjuntos de datos |
| Random Forests | ensemble.RandomForestsClassifier |
Bagging de árboles de decisión |
| Extremelly Randomized Trees | ensemble.ExtraTreesClassifier |
Lanza árboles randomizados |
| ADABoost | ensemble.AdaBoostClassifier |
Adaptative Boosting |
| Gradient Boosted Trees | ensemble.GradientBoostingClassifier |
Boosting de árboles |
Vamos a explorar algunos:
4.A. Bagging¶
Vamos a probar a hacer un bagging de regresiones logísticas sobre nuestro dataset sobre la alopecia, utilizando el BaggingClassifier:
from sklearn.ensemble import BaggingClassifier
# Generemos nuestra regresión logística:
lr = LogisticRegression(penalty="l1", # Penalización L1 o LASSO
solver="liblinear")
# Y ahora generamos el bagging:
bagging_lr = BaggingClassifier(lr, # Modelos
n_estimators=50, # Número de modelos en el bagging
max_samples=0.3, # Proporción máxima del dataset dada a cada modelo para entrenar
max_features=1, # Número máximo de features dado a cada modelo para entrenar
bootstrap=True, # Con reemplazamiento en la parte del dataset cogida por modelo
bootstrap_features=True, # Con reemplazamiento en las features
n_jobs=4 # Paralelizar el trabajo en 4 núcleos
)
Al igual que con un modelo normal y corriente, lo entrenamos con .fit():
bagging_lr.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
Y comprobamos el rendimiento en el conjunto de testing:
print("Bien clasificados en testing: %.16f" % bagging_lr.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"]))
figura = plt.figure(figsize=(16,7))
figura.suptitle("Bagging de LR (50 modelos)", fontsize=14, fontweight="bold")
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(bagging_lr,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(bagging_lr,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
Vemos que lo hace muy bien en el conjunto de testing. Aunque no se aprecia, el resultado es no lineal. Parece que el modelo generaliza muy bien.
4.B. Random Forests¶
Random Forests es uno de los modelos más famosos para clasificación. No deja de ser un bagging de árboles de decisión, con una pequeña modificación (básicamente que cada árbol es entrenado con subconjuntos no solo de las observaciones, sino también de las features). Pero un bagging de árboles al fin y al cabo. Suele dar muy buenos resultados.
Probemos entonces el RandomForestClassifier:
from sklearn.ensemble import RandomForestClassifier
rf = RandomForestClassifier(n_estimators=20, # 20 árboles
criterion="entropy", # Igual que en el DecisionTreeClassifier
max_depth=3, # Igual que en el DecisionTreeClassifier
min_samples_split=10,
min_samples_leaf=5,
bootstrap=True)
rf.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
print("Bien clasificados en testing: %.16f" % rf.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"]))
figura = plt.figure(figsize=(16,7))
figura.suptitle("Random Forests (20 árboles)", fontsize=14, fontweight="bold")
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(rf,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(rf,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
4.C. AdaBoost¶
AdaBoost (o Adaptive Boosting) es el algoritmo de boosting más popular y, en muchas ocasiones, el que mejores resultados da.
AdaBoost va pasando clasificadores (o regresores) sobre formas modificadas de los datos, de forma repetida. Por cada boosting iteration (pasada de un nuevo clasificador), se asignan pesos o importancias a cada punto, de forma que los puntos mal clasificados por el modelo anterior cobran más importancia; mientras que los puntos bien clasificados la van perdiendo. Esta técnica adaptativa de boosting es muy popular y bastante efectiva en general.
Utilicemos sklearn.ensemble.AdaBoostClassifier para generar un adaptive boosting de regresiones logísticas:
from sklearn.ensemble import AdaBoostClassifier
lr = LogisticRegression(penalty="l2",
solver="lbfgs")
ab = AdaBoostClassifier(base_estimator=lr,
n_estimators=20, # Número de clasificadores
learning_rate=0.4 # La "contribución" de cada clasificador
)
ab.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
print("Bien clasificados en testing: %.16f" % ab.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"]))
figura = plt.figure(figsize=(16,7))
figura.suptitle("AdaBoost (Logistic Regression)", fontsize=14, fontweight="bold")
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(ab,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(ab,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
4.D. Gradient Boosting¶
La técnica de Gradient Boosting es, cuanto menos, curiosa: al igual que el resto de técnicas de boosting, se lanza primero un clasificador (en el caso de Gradient Boosting suelen ser árboles de decisión). Luego, se formula el árbol generado como una función diferenciable; y optimizable (mediante la computación de gradientes e ir descendiendo el error generado) a través del árbol generado para la siguiente iteración. Los señores de XGBoost (biblioteca que mencionaremos al final) lo explican mejor en su blog.
De hecho, XGBoost es una implementación de Gradient Boosting (no incluida en scikit-learn, sino como una biblioteca de Python aparte) escalable y que incorpora una buena regularización (luego hablaremos de las regularizaciones). Pero aquí vamos a utilizar la implementación de Gradient Boosting de scikit-learn.
Lancemos un GradientBoostingClassifier:
from sklearn.ensemble import GradientBoostingClassifier
gbt = GradientBoostingClassifier(learning_rate=0.05,
n_estimators=40,
max_depth=3, # Son árboles al fin y al cabo
min_samples_split=5,
min_samples_leaf=5
)
gbt.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
print("Bien clasificados en testing: %.16f" % gbt.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_testing["label"]))
figura = plt.figure(figsize=(16,7))
figura.suptitle("Gradient Boosted Trees (30 Trees)", fontsize=14, fontweight="bold")
# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(gbt,
dataset_clasificacion_training[["hormona_a", "hormona_b"]],
dataset_clasificacion_training["label"],
title="Training",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes1)
# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(gbt,
dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
dataset_clasificacion_testing["label"],
title="Testing",
xlabel="Concentración de hormona a",
ylabel="Concentración de hormona b",
fig=axes2)
Los ensembles de modelos se utilizan mucho en clasificación; pero también se pueden utilizar en problemas de regresión: apartado que veremos después de las métricas de clasificación.
5. Métricas de clasificación¶
Hasta ahora hemos utilizado una métrica para medir el rendimiento de nuestros clasificadores. Al aplicar el método .score() en un clasificador entrenado, nos devuelve un número entre 0 y 1. Este número nos dice el porcentaje de observaciones (o puntos) bien clasificadas; métrica conocida como accuracy.
Esta métrica está muy bien en algunos casos, pero no siempre. Supongamos un dataset como el siguiente (vamos a crear solo las labels, por simplificar):
np.concatenate((np.zeros(90),np.ones(10)))
Donde tenemos 90 observaciones que son de la clase cero, y solo 10 que son de la clase 1. Pues bien, supongamos que creamos un "modelo" como el siguiente:
Nuestro modelo predice siempre 0, independientemente de la observación.
Con este modelo tan estúpido conseguiríamos un 90% de puntos correctamente clasificados. No obstante, eso no quiere decir que sea un buen modelo (evidentemente no lo es). Más importante que clasificar bien un X% de los puntos es que nuestro modelo sea capaz de discriminar, y diferenciar lo que debería ser un uno de lo que debería ser un cero.
Por esa razón la métrica de porcentaje de puntos bien clasificados no es siempre ideal.
5.A. Confusion matrix¶
Afortunadamente, existen otras métricas que nos permiten ver el rendimiento de un clasificador.
Una de las primeras cosas que podemos hacer tras entrenar un clasificador es construir lo que se conoce como confusion matrix o matriz de confusión. La matriz de confusión no es más que una tablita que se utiliza normalmente para ver qué tal lo ha hecho un clasificador en problemas de clasificación binaria, es decir, problemas de clasificar dos labels (aunque se puede extender a problemas de clasificar más de dos clases). Esta tablita tiene cuatro cuadrantes en el caso de la clasificación binaria:
| Predicho 1 | Predicho 0 | |
|---|---|---|
| Valor real 1 | True Positives | False Negatives |
| Valor real 0 | False Positives | True Negatives |
- El cuadrante superior izquierdo muestra los true positives, es decir, el número de unos bien clasificados.
- El cuadrante inferior izquierdo muestra los false positives: el número de ceros clasificados como uno.
- El cuadrante superior derecho muestra los false negatives: el número de unos clasificados como cero.
- El cuadrante inferior derecho muestra los true negatives: el número de ceros bien clasificados.
Esta matriz nos permite, de un rápido vistazo, determinar si estamos "sobreclasificando" de alguna de las dos clases (por ejemplo: muchísimos ceros clasificados como unos, pero pocos unos clasificados como ceros).
sklearn.metrics incluye la función confusion_matrix que nos permite construir la matriz de confusión. Toma como primer argumento las labels reales que vemos en las observaciones, y como segundo argumento las labels predichas por un modelo (las cuales podemos conseguir con modelo.predict()).
Dicho esto, vamos a construir la matriz de confusión del último clasificador que hemos entrenado, el gradient boosted trees que hemos llamado gbt. Utilizaremos los datos del conjunto de testing para ello:
from sklearn.metrics import confusion_matrix
matriz_confusion = confusion_matrix(y_true=dataset_clasificacion_testing["label"],
y_pred=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
)
# Sklearn nos saca en la primera columna los
# predichos como cero, en vez de los predichos
# como uno, así que cambiamos el orden de las
# columnas y filas para que sea como el de nuestra
# definición previa:
matriz_confusion = np.array([(matriz_confusion[1]),
matriz_confusion[0]])[:,::-1]
matriz_confusion
Acorde a los cuatro cuadrantes definidos antes, nuestro gbt en el conjunto de testing:
- Ha predicho bien 55 unos -> true positives
- Ha predicho 2 ceros como uno (2 hombres no calvos como calvos) -> false positives
- Ha predicho 6 unos como cero (6 hombres calvos como no calvos) -> false negatives
- Ha predicho bien 57 ceros -> true negatives
En el gráfico que hicimos cuando entrenamos el modelo resulta fácil de verificar.
Podemos ver que no existe ningún "fallo gordo" con el modelo. El problema vendría por ejemplo si las cantidades de false positives y de false negatives fueran muy dispares.
5.B. Precision, recall y F1-score¶
La matriz de confusión es fácil de interpretar, pero para nosotros; no para los ordenadores. Ellos prefieren un número, o un par de ellos más fáciles de comparar entre modelos.
Para eso existen precisamente otras métricas como son la precision, el recall y el F1-score:
5.B.a. Precision¶
La precision puede entenderse como la habilidad que tiene un clasificador para no predecir los ceros como unos. Su cálculo es muy sencillo:
$$precision=\frac{\sum true\,positives}{\sum (true\,positives + false\,positives)}$$O dicho de otro modo:
$$precision=\frac{\sum Los\,adecuadamente\,predichos\,como\,1}{\sum Todos\,los\,predichos\,como\,1}$$Podemos computar fácilmente la precision con sklearn.metrics.precision_score. Se utiliza exactamente igual que confusion_matrix, pero en este caso devuelve un solo número en vez de una matriz:
from sklearn.metrics import precision_score
precision_gbt = precision_score(y_true=dataset_clasificacion_testing["label"],
y_pred=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
)
print("Precision: " + str(precision_gbt))
print("Cuando nuestro clasificador predice una observación como uno, acierta en el "+ str(precision_gbt*100)
+ "% de los casos (en el conjunto de testing).")
5.B.b. Recall¶
El recall es la habilidad que tiene un clasificador para predecir todos los unos reales. Su cálculo es:
$$recall=\frac{\sum true\,positives}{\sum (true\,positives + false\,negatives)}$$O bien:
$$recall=\frac{\sum Los\,adecuadamente\,predichos\,como\,1}{\sum Todos\,los\,que\,son\,1\,en\,realidad}$$Computamos el recall con sklearn.metrics.recall_score:
from sklearn.metrics import recall_score
recall_gbt = recall_score(y_true=dataset_clasificacion_testing["label"],
y_pred=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
)
print("Recall: " + str(recall_gbt))
print("Nuestro clasificador ha conseguido predecir bien el " + str(recall_gbt*100)
+ "% de los unos (en el cojunto de testing).")
5.B.c. F1-score¶
El F1-score esla métrica más utilizada dentro de las que son llamadas F-measures (valores-F en castellano). Básicamente combina en una sola métrica la precision y el recall a través de la siguiente fórmula:
$$F_1\mbox{-}score=2 \times\frac{precision \times recall}{precision + recall}$$El F1-score da la misma importancia a la precision y al recall. Otras F-measures pueden utilizarse para dar más importancia a una métrica que a la otra. Qué F-measure utilizar depende mucho del negocio y del problema a resolver.
Utilicemos sklearn.metrics.f1_score para computarlo:
from sklearn.metrics import f1_score
f1_score_gbt = f1_score(y_true=dataset_clasificacion_testing["label"],
y_pred=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
)
print("F1-score: " + str(f1_score_gbt))
Como regla general: cuanto más cercano sea el F1-score a 1, mejor; porque querrá decir que la precision y el recall también son cercanos a 1.
5.C. Curva ROC¶
La curva ROC (Receiver Operating Characteristic) es un gráfico muy útil para visualizar el rendimiento de un clasificador binario. Este gráfico tiene los siguientes ejes:
- El eje x es el False Positive Rate (FPR).
- El eje y es el True Positive Rate (TPR).
Donde el False Positive Rate es:
$$False\,Positive\,Rate=\frac{\sum false\,positives}{\sum (false\,positives+true\,negatives)}$$O:
$$False\,Positive\,Rate=\frac{\sum Los\,0\,mal\,predichos\,como\,1}{\sum Todos\,los\,que\,son\,0\,en\,realidad}$$Y el True Positive Rate es lo mismo que el recall:
$$True\,Positive\,Rate\,o\,recall=\frac{\sum true\,positives}{\sum (true\,positives + false\,negatives)}$$Esta curva tiene dos usos principales:
- En el caso de clasificadores que generan un scoring, y que luego determinamos que son ceros las observaciones que puntúan por debajo de cierto valor umbral; y unos los que puntúan por encima: la curva ROC nos permite visualizar cuál es el mejor umbral que podemos poner. Por ejemplo: la regresión logística nos devuelve valores entre 0 y 1 para cada observación, y por defecto determinamos que los que den < 0.5 serán predichos como cero, y que los que den >= 0.5 serán predichos como uno. Pues bien, la curva ROC nos permite ver si existe otro umbral (0.7, 0.4...) que sea mejor para nuestro problema, en función de lo que busquemos. En estos casos, la curva ROC estará formada por los puntos de (FPR, TPR) para distintos umbrales.
- En el caso de clasificadores que no generen scoring la "curva ROC" será solo una línea que une el (0,0), el punto (FPR, TPR) de nuestro modelo, y el (1,1). Medir el área bajo esta curva (llamada Área Under ROC o AUC) es una métrica de rendimiento que nos permite comparar clasificadores (también se puede usar así para clasificadores que sí generan scoring).
Aquí nos centraremos en el segundo uso. Puedes aprender más sobre el primero con este vídeo.
Lo ideal es que un modelo tenga un False Positive Rate muy bajo (lo más cercano a $0$ posible), y un True Positive Rate muy elevado (lo más cercano a $1$ posible). De esta forma, el área bajo la curva ROC (AUC) será lo más elevada posible. Veámoslo con unos ejemplos.
Supongamos que hacemos un clasificador, y obtenemos que su FPR es $0.3$; y su TPR es $0.4$. Con este punto, el (0,0) y el (1,1) podemos dibujar la curva ROC:
coordenadas_fpr = [0, 0.3, 1]
coordenadas_tpr = [0, 0.4, 1]
plt.figure(figsize=(8,6))
# Primero pintamos el random guess:
plt.plot([0,1],[0,1], "g--")
# Y ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")
# Coloreamos el área bajo la curva
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.4, color="yellow")
plt.title("Curva ROC (modelo mediocre)")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
pass
La línea verde quebrada muestra lo que se llama como random guess o línea de no-discriminación, es decir: lo que haría un modelo que predice 0 o 1 "a boleo". Este modelo acertaría aproximadamente el 50% de los casos, dando por hecho que tenemos el mismo número de ceros que de unos. Se suele representar para comparar con la curva ROC de nuestro modelo.
Y la línea roja es la curva ROC de este modelo. Este modelo es poco mejor que uno que adivinara "a boleo": lo ideal es que la curva crezca lo más posible al principio, "despegándose" lo más posible de la línea de random guess por arriba. El área bajo esta curva ROC (coloreada de amarillo) es cercana a $0.5$. Teniendo en cuenta que el área total del gráfico es $1$, es poco más de la mitad: lo cual no es ideal. Fijándonos en el punto (0.3,0.4) de nuestro gráfico, podemos interpretar esta curva ROC como:
Para detectar y acertar el 40% de los unos, me "trago" o clasifico mal un 30% de los ceros de mi dataset.
Un clasificador perfecto tendría un FPR de $0$ y un TPR de $1$; y el área bajo su curva ROC sería $1$ (es decir, todo el gráfico).
Otro ejemplo: ahora supongamos que entrenamos otro modelo, y éste tiene un FPR también de 0.3; pero que por el contrario su TPR es más alto, siendo 0.8:
coordenadas_fpr = [0, 0.3, 1]
coordenadas_tpr = [0, 0.8, 1]
plt.figure(figsize=(8,6))
# Primero pintamos el random guess:
plt.plot([0,1],[0,1], "g--")
# Y ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")
# Coloreamos el área bajo la curva
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.4, color="yellow")
plt.title("Curva ROC (modelo mejor)")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
pass
Este está mejor: nuestra curva ROC ahora está bastante "más alta" que el random guess, y el AUC es bastante mejor. Sigue sin estar cerca de ser $1$, pero definitivamente es mayor; y bastante mayor que $0.5$.
Para detectar y acertar el 80% de los unos, me "trago" o clasifico mal un 30% de los ceros de mi dataset.
Y por último, supongamos un modelo muy malo, con FPR de 0.7 y TPR de 0.3:
coordenadas_fpr = [0, 0.7, 1]
coordenadas_tpr = [0, 0.3, 1]
plt.figure(figsize=(8,6))
# Primero pintamos el random guess:
plt.plot([0,1],[0,1], "g--")
# Y ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")
# Coloreamos el área bajo la curva
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.4, color="yellow")
plt.title("Curva ROC (modelo penoso)")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
pass
Este modelo es claramente pésimo: hace peor trabajo que el random guess, es decir: es mejor predecir "a boleo" que utilizar este modelo. Su AUC es bastante menor a $0.5$.
Para detectar y acertar el 30% de los unos, me "trago" o clasifico mal un 70% de los ceros de mi dataset.
Para usar un modelo como este, mejor nos dedicamos al arte abstracto.
Ahora, probemos a dibujar la curva ROC de nuestro modelo gbt.
sklearn.metrics.roc_curve nos genera las coordenadas para nuestros modelos de forma muy similar a como las hemos metido en los ejemplos que acabamos de inventarnos. La función nos devuelve una tupla con tres arrays de Numpy:
- El primer array nos da las coordenadas del eje x, es decir: del FPR.
- El segundo array nos da las coordenadas del eje y, es decir: el TPR.
- El tercer array nos da los thresholds o umbrales, en el caso de que quisiéramos utilizar la curva ROC para elegir el mejor umbral para un modelo de scoring (acorde al que pusimos como primero de los dos posibles usos para las curvas ROC). Así que aquí vamos a ignorarla: no tiene sentido utilizarla sin modelo de scoring.
from sklearn.metrics import roc_curve
roc = roc_curve(y_true=dataset_clasificacion_testing["label"],
y_score=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]]),
pos_label=1.0 # Le decimos a roc_curve que los casos positivos son los unos.
)
roc
Ahora cogemos las dos primeras arrays del resultado, y dibujamos la curva ROC de nuestro modelo gbt:
# Cogemos las coordenadas del resultado:
coordenadas_fpr = roc[0]
coordenadas_tpr = roc[1]
plt.figure(figsize=(8,6))
# Primero pintamos el random guess:
plt.plot([0,1],[0,1], "g--")
# Y ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")
# Coloreamos el área bajo la curva
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.4, color="yellow")
plt.title("Curva ROC modelo gbt")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
pass
Creo que queda claro que nuestro modelo gbt es muy bueno: la curva ROC se "despega" muchísimo del random guess por arriba; y crece muchísimo al principio. Por eso mismo, el AUC es casi toda la posible; y debe ser muy cercana a $1.0$.
De hecho, sklearn.metrics.roc_auc_score nos permite calcular fácilmente el AUC. Se utiliza como cualquier otra métrica de clasificación:
from sklearn.metrics import roc_auc_score
auc_gbt = roc_auc_score(y_true=dataset_clasificacion_testing["label"],
y_score=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
)
print("AUC: " + str(auc_gbt))
Efectivamente: el área bajo la curva ROC es muy cercana al máximo, que es $1$. Concluimos: AUC es una métrica sencilla de utilizar para evaluar la efectividad de nuestros clasificadores (al igual que la precision, recall y F1-score); y dibujar la curva ROC también es una forma vistosa de comparar clasificadores.
Existen muchas más métricas para medir el rendimiento de clasificadores, si bien las mencionadas son las más famosas y utilizadas. Otros ejemplos pueden ser la curva Lift (muy utilizada en marketing) o el LogLoss (disponible en sklearn).
Y hay todavía más, disponibles en sklearn.metrics: tanto para clasificación como para regresión y otros.
Ya está bien de clasificación. Pasemos al otro tipo principal de problema supervisado: la regresión.
6. Regresión¶
A diferencia de la clasificación (donde predecimos labels o clases), en regresión predecimos un número; también en función de una serie de features de cada observación. Un ejemplo podría ser el visto al principio de este documento: intentar predecir el precio de las casas en Boston a partir de una serie de variables.
Otro ejemplo podría ser predecir la renta de una persona a partir de features, que es el ejemplo que vamos a hacer. El caso es que en regresión la label (el valor a predecir) se le suele llamar target, y debe ser una variable continua (numérica, y a ser posible con decimales). Todos los conceptos de sobreajuste, generalización y similares aplican también a los problemas de regresión.
Sin más dilación, comencemos por generar un dataset sintético para regresión. sklearn contiene la función sklearn.datasets.make_regression, pero los datasets que genera son demasiado sencillos: suele ser muy fácil ajustar el problema de regresión. Así que vamos a usar álgebra con Numpy para generar un dataset más interesante para regresión.
En los problemas de clasificación, como hemos visto hasta ahora, se puede visualizar bien "3 dimensiones en 2D". Quiero decir: teníamos dos features y la label, lo cual en realidad suman 3 dimensiones distintas. Lo visualizábamos dibujando en 2D las features, y poniendo la label como color, a modo de la tercera dimensión.
En regresión suele ser más complicado visualizarlo así, así que vamos a generar solo una feature que irá en el eje x, el target que irá en el eje y, y el color de los puntos será uniforme e indiferente. Vamos a ello: pondremos en el eje x los años trabajando de distintos individuos, y en el eje y el target, que será el salario bruto anual:
# Comencemos por generar la coordenada x:
np.random.seed(12)
x = np.random.normal(loc=9, scale=16, size=3200)
# Vamos a quitar todos los valores que sean menores de 0
# (porque nadie puede llevar -1 años trabajados):
x = x[x>=0]
# Nos han quedado estas observaciones:
print("Número de observaciones: " + str(len(x)))
# Y ahora la coordenada y, es decir, el target,
# a modo de función de x. Y vamos a añadir
# algo de ruido metiendo algo de random:
y = 8000*x**(1.0/3.0) + 0.8*x**3 + 21000 + np.random.normal(scale=5000, size=len(x))
# A DataFrame:
dataset_regresion = pd.DataFrame({"años trabajados": x, "salario bruto": y})
# Primeras cinco filas de nuestro dataset de regresión:
dataset_regresion[:5]
dataset_regresion.plot.scatter(x="años trabajados",y="salario bruto",figsize=(8,6))
pass
Dividimos en train y test:
dataset_regresion_spliteado = train_test_split(dataset_regresion,
train_size=0.8, # 80% training
test_size=0.2, # 20% testing
random_state=24 # seed para el generador de num. al.
)
dataset_regresion_training = dataset_regresion_spliteado[0]
dataset_regresion_testing = dataset_regresion_spliteado[1]
Y los dibujamos:
plt.figure(figsize=(16,7))
# Training
axes1=plt.subplot(1,2,1, title="Training")
dataset_regresion_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)
# Testing
axes2=plt.subplot(1,2,2, title="Testing")
dataset_regresion_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
pass
No es súper realista, pero nos servirá para nuestros propósitos.
Sklearn incluye muchos modelos regresores. De hecho, prácticamente todos los modelos clasificadores tienen su variante regresora. Estos son algunos de ellos:
| Algoritmo | Módulo | Comentario |
|---|---|---|
| Mínimos Cuadrados | linear_model.LinearRegression |
Regresión lineal con Least Squares |
| Ridge Regression | linear_model.Ridge |
Least Squares con regularización $L2$ |
| Lasso Regression | linear_model.Lasso |
Least Squares con regularización $L1$ |
| Elastic Net | linear_model.ElasticNet |
Least Squares con regularización $L1$ y $L2$ |
| KNN Regression | neighbors.KNeighborsRegressor |
KNN adaptado a regresión (filtro de media) |
| Radius NN Regression | neighbors.RadiusNeighborsRegressor |
NN teniendo en cuenta el radio y no el número de vecinos de la muestra |
| Regression Trees | tree.DecisionTreeRegressor |
Árboles de decisión para regresión |
| Linear SVM Regresión | svm.LinearSVR |
Implementación de SVM para regresión. Escala mejor para un gran número de muestras. Emplea la biblioteca liblinear |
| SVM Regression | svm.SVR |
Implementación de SVM para regresión. Admite kernels. Emplea la biblioteca libsvm |
| Isotonic Regression | isotonic.IsotonicRegression |
Regresión diseñada para datos con valor creciente |
| Procesos Gaussianos | gaussian_process.GaussianProcess |
Método de regresión con salida probabilística |
| Regresion Polinomica | preprocessing.PolynomialFeatures |
Expansión de los datos para regresión no lineal |
Comencemos a hacer regresiones:
6.A. Regresión lineal simple (mínimos cuadrados ordinarios)¶
En el ejercicio de Numpy aprendimos a hacer con álgebra lineal el ajuste por mínimos cuadrados. Sklearn nos permite hacerlo de forma sencilla con sklearn.linear_model.LinearRegression:
from sklearn.linear_model import LinearRegression
ols = LinearRegression(fit_intercept=True, # Añadir intercept, bias u ordenada en el origen
normalize=False, # Normalizar las features antes de entrenar. Es el valor por
) # defecto, pero lo pongo para que lo veas.
Exactamente igual que con los clasificadores, llamamos al método .fit() para entrenar la regresión:
ols.fit(X=dataset_regresion_training[["años trabajados"]], # Aunque sea una sola feature,
y=dataset_regresion_training["salario bruto"]) # sklearn quiere que la seleccionemos
# como una lista de columnas, con solo
# un elemento
El problema de regresión lineal planteado es el siguiente:
$$salario = \beta_{0} + \beta_{1}años$$Ya hemos entrenado el modelo, así que tenemos los parámetros $\beta$ estimados:
print("beta0 o intercept: " + str(ols.intercept_))
print("beta1: " + str(ols.coef_[0]))
El rendimiento en regresión no se puede medir en acertar/fallar el label, puesto que se trata de números continuos. Así, podemos medir cuánto nos hemos acercado con nuestra predicción al número real. Sklearn tiene en sus regresores (una vez entrenados) el método .score(), exactamente igual que en los clasificadores.
No obstante, cuando aplicamos .score(), sklearn calcula el $R^2$ (coeficiente de determinación), el cual NO recomiendo como medida del rendimiento de nuestra regresión. Vamos a probarlo de todos modos:
print("R cuadrado : " + str(ols.score(X=dataset_regresion_testing[["años trabajados"]],
y=dataset_regresion_testing["salario bruto"])))
En problemas de regresión, recomiendo utilizar otro tipo de métricas, como puede ser el error cuadrático medio (mean squared error o $MSE$), o bien la raíz de este valor (root mean squared error o $RMSE$). Afortunadamente, sklearn tiene varias métricas para medir el rendimiento de regresores en sklearn.metrics.
Vamos a calcular el $MSE$ y $RMSE$ utilizando sklearn.metrics.mean_squared_error:
# mean_squared_error toma como primer argumento los valores
# reales del target, y como segundo los valores que predecimos:
from sklearn.metrics import mean_squared_error
# Medimos el MSE en el conjunto de testing.
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
y_pred=ols.predict(dataset_regresion_testing[["años trabajados"]]))
# El RMSE es la raíz de este valor. Podemos utilizar
# np.sqrt o math.sqrt (de la biblioteca estándar)
# para calcularlo:
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " + str(rmse))
Y este es nuestro $E_{out}$ para nuestra regresión lineal. El $RMSE$ se puede entender como: en términos medios, fallamos las predicciones por unos 10.600€ de salario (aprox), ya sea de más o de menos, sobre el valor real.
También podemos dibujar nuestra línea de regresión:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Regresión lineal simple (mínimos cuadrados ordinarios)", fontsize=14, fontweight="bold")
# Training
axes1=plt.subplot(1,2,1, title="Training")
dataset_regresion_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)
# Dibujamos la recta de regresión:
plt.plot(dataset_regresion_training["años trabajados"],
ols.predict(dataset_regresion_training[["años trabajados"]]),
c="red")
# Testing
axes2=plt.subplot(1,2,2, title="Testing")
dataset_regresion_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
# Dibujamos la recta de regresión:
plt.plot(dataset_regresion_testing["años trabajados"],
ols.predict(dataset_regresion_testing[["años trabajados"]]),
c="red")
pass
El ajuste no es ideal. La regresión lineal simple no puede hacer mucho más...
Teniendo en cuenta la "forma" de nuestro dataset, ningún modelo lineal va a poder darnos resultados mucho mejores. Dicho esto, vamos a mencionar brevemente tres de las regularizaciones más famosas para la regresión lineal:
6.B. Regularización: Ridge, LASSO y Elastic Net¶
En clasificación vimos que algunos modelos toman como argumento algunos parámetros como "l1" o "l2". Esto son lo que se llaman regularizaciones. Otros eran simplemente argumentos que especificaban las características del modelo (como puede ser el número de vecinos a contemplar en KNN, o la profundidad del árbol en los árboles de decisión).
La regularización es un proceso que intenta alterar ligeramente la formulación matemática de un modelo de machine learning, con la intención de prevenir el sobreajuste. Realmente, consiste en añadir hiperparámetros a nuestro modelo que son elegidos por nosotros, que incluyen ciertas penalizaciones. Dichas penalizaciones hacen que cuando el modelo se entrene, lo haga acorde a reducir el error (como siempre). Pero ahora este error ha sido alterado, incluyendo dichas penalizaciones; de forma que el modelo intentará conseguir un balance entre ajustar bien, y no ser penalizado en exceso.
Diseñar una regularización es un proceso algo complicado: requiere saber qué se quiere penalizar, cómo añadirlo a la formulación matemática del modelo, y medir los resultados. Dicho esto, existen tantas posibles regularizaciones como las que se nos pueda ocurrir.
Tanto para la regresión lineal como para la regresión logística (y también aplicable a otros varios modelos) existen una serie de regularizaciones que se conoce que funcionan bastante bien y están bastante consolidadas. Las más famosas y utilizadas son estas dos:
- Regularización o penalización $L2$: no entraré en terreno matemático. A efectos prácticos, la regularización $L2$ consigue que los parámetros estimados por nuestro modelo (las $\beta$ según la nomenclatura que hemos ido siguiendo) no tengan un valor absoluto demasiado grande. Es decir: introducir la penalización $L2$ en una regresión lineal hará que el modelo "sea forzado" a estimar parámetros $\beta$ pequeños (0.5, 3.1, -0.6) antes que parámetros $\beta$ grandes (3.320, 4.930, -14.291). Debido a esto, esta regularización se conoce también como parameter shrinking (contracción de parámetros estimados; puesto que los hace chiquititos) . El nombre "científico" es regularización de Tikhonov.
- Regularización o penalización $L1$: es más curiosa aún. Consigue que algunos de los parámetros $\beta$ estimados sean 0, es decir, que no sean útiles ni afecten a la predicción. Así que en el fondo, lo que hace es forzar la selección de features o variables dentro del modelo. Que yo sepa, no tiene nombre "científico".
Como estas dos regularizaciones son tan famosas, las regresiones lineales que las aplican tienen nombres específicos (aunque repito que se pueden utilizar en modelos más allá de las regresiones lineales):
- Regresión Ridge: es la regresión lineal que aplica la regularización $L2$.
- Regresión LASSO (o LASSO a secas): es la regresión lineal que aplica la regularización $L1$.
- Elastic Net: la regresión lineal que aplica tanto $L1$ como $L2$, y nos permite fijar qué proporción de cada una de ellas aplica.
A cada una de estas regularizaciones les acompaña por tanto un hiperparámetro que fijamos nosotros, y que dicta la "intensidad" de estas regularizaciones/penalizaciones.
Mucha teoría. Vamos a hacer un ejemplo utilizando la regresión Ridge, disponible en sklearn en sklearn.linear_model.Ridge:
from sklearn.linear_model import Ridge
# Formulamos la Ridge:
ridge = Ridge(alpha=1, # Hiperparámetro de la regularización, que marca la intensidad del shrinking
fit_intercept=True)
# La entrenamos:
ridge.fit(X=dataset_regresion_training[["años trabajados"]],
y=dataset_regresion_training["salario bruto"])
Vamos a ver si los parámetros $\beta$ estimados efectivamente se han hecho más pequeños, comparados con los de la regresión lineal simple por mínimos cuadrados que hicimos antes. Prestamos atención a $\beta_1$ (y no al $\beta_0$ o intercept):
print("beta0 o intercept: " + str(ridge.intercept_))
print("beta1: " + str(ridge.coef_[0]))
Algo lo ha reducido, pero no mucho. Si aumentamos el hiperparámetro alpha, es decir, la intensidad de la regularización:
# Formulamos la Ridge:
ridge2 = Ridge(alpha=100000, # No hagas esto nunca; es demsiado alto...
fit_intercept=True)
# La entrenamos:
ridge2.fit(X=dataset_regresion_training[["años trabajados"]],
y=dataset_regresion_training["salario bruto"])
print("beta0 o intercept: " + str(ridge2.intercept_))
print("beta1: " + str(ridge2.coef_[0]))
Efectivamente, lo va reduciendo. En el caso de que tuviéramos más features (y por tanto más estimadores), el efecto sería mucho más notable. Cuando tenemos solo una variable, la regularización no surte tanto efecto; de forma que lo que hemos hecho aquí es poco útil, y solo tiene propósitos educativos.
Siéntete libre de probar también con LASSO y Elastic Net para ver cómo cambian los parámetros estimados, a ser posible con un problema multivariante (es decir, con más de una feature) para ver mejor los efectos de las regularizaciones.
La pregunta que queda es: ¿cómo elegir una regularización, así como la intensidad de la misma? La respuesta es la misma que cuando se trata de seleccionar modelos (al fin y al cabo, sigue siendo seleccionar el modelo, o su variante regularizada, que mejor funciona). Para ello, veremos luego el proceso de validación de modelos y su selección.
Ahora, continuamos con las regresiones. Vamos a probar la polinómica.
6.C. Regresión polinómica: transformación no lineal de nuestro dataset¶
Las regularizaciones sobre un modelo lineal siguen dándonos un modelo lineal. Vamos a probar a hacer regresiones no lineales utilizando polinomios. sklearn.preprocessing incluye las PolynomialFeatures, que básicamente coge nuestras features y las transforma en polinomios del grado que le pidamos.
Normalmente esta función es más útil cuando tenemos más de una feature. Supongamos que tenemos dos (años trabajados y altura de la persona):
$$salario = \beta_{0} + \beta_{1}años + \beta_{2}altura$$Aplicando PolynomialFeatures de grado 2, nuestro problema se transformaría en:
$$salario = \beta_{0} + \beta_{1}años + \beta_{2}altura + \beta_{3}años^2 + \beta_{4}altura^2 + \beta_{5}(años\times altura)$$
Lo cual nos podría dar un mejor ajuste, ya que contempla mayor nivel de interacción entre las features.
En nuestro caso donde solo tenemos una feature, la expansión polinómica de grado 2 generará:
$$salario = \beta_{0} + \beta_{1}años + \beta_{2}años^2$$Vamos a probarlo:
from sklearn.preprocessing import PolynomialFeatures
# Ahora, creamos una instancia de PolynomialFeatures
# de grado 2:
poly2 = PolynomialFeatures(degree=2,
include_bias=False) # Ya añadiremos el bias/intercept cuando entrenemos modelo
# Y ahora transformamos nuestro dataset, utilizando
# un método para los transformadores de features
# llamado fit_transform():
features_poly2 = poly2.fit_transform(dataset_regresion_training[["años trabajados"]])
A partir de lo que solo era una feature (años trabajados), se han creado dos:
features_poly2
Vamos a meterlas en un DataFrame nuevo:
dataset_poly2_training = pd.DataFrame(features_poly2)
# Nombramos columnas:
dataset_poly2_training.columns = ["años trabajados", "años trabajados ^2"]
# Y añadimos los targets del conjunto training original, inalterados.
# OJO: queremos que pandas "pase" de los índices porque nos desordenaría
# la columna de los targets con respecto al resto, así que
# quitamos los índices pasando la columna a array de Numpy:
dataset_poly2_training["salario bruto"] = np.array(dataset_regresion_training["salario bruto"])
Así era el conjunto de training antes de las PolynomialFeatures:
dataset_regresion_training[:7]
Y así es después:
dataset_poly2_training[:7]
Muy bien. Ahora entrenemos la regresión. El trabajo de convertir la regresión lineal en regresión polinómica ya está hecho: ya hemos aplicado la transformación no lineal a nuestra feature:
Aplicar una regresión (o modelo) lineal en un espacio dimensional (o dataset) resultado de la transformación no lineal de las features, da como resultado un modelo no lineal.
Entrenamos (regresión lineal simple de mínimos cuadrados ordinarios):
lr_2 = LinearRegression(fit_intercept=True)
lr_2.fit(X=dataset_poly2_training[["años trabajados", "años trabajados ^2"]],
y=dataset_poly2_training["salario bruto"])
Hecho. Ahora, no obstante, tenemos también que aplicar la transformación al conjunto de testing, ya que para predecir, tenemos igualmente que realizar la transformación polinómica a las features:
# Hacemos fit_transform del PolynomialFeatures
# y metemos en DataFrame nuevo:
dataset_poly2_testing = pd.DataFrame(poly2.fit_transform(dataset_regresion_testing[["años trabajados"]]))
# Renombramos columnas:
dataset_poly2_testing.columns = ["años trabajados", "años trabajados ^2"]
# Y añadimos los targets:
dataset_poly2_testing["salario bruto"] = np.array(dataset_regresion_testing["salario bruto"])
Conjunto de testing antes de la transformación:
dataset_regresion_testing[:7]
Y después:
dataset_poly2_testing[:7]
Medimos el $RMSE$ en el conjunto de testing:
rmse_poly2 = np.sqrt(mean_squared_error(y_true=dataset_poly2_testing["salario bruto"],
y_pred=lr_2.predict(dataset_poly2_testing[["años trabajados","años trabajados ^2"]])))
print("RMSE: " + str(rmse_poly2))
Parece que hemos conseguido dejar el $E_{out}$ a la mitad con nuestra regresión polinómica (comparando con la regresión lineal simple). Dibujemos:
# Cuidado: para dibujar curvas, es importante
# que los datos en el eje x estén ordenados de menor a mayor
# (y por tanto, re-ordenar también los del)
# eje y acorde. En este caso, como las coordenadas
# de nuestras curvas se obtienen con lr_2.predict(),
# con ordenar los datos del eje x (años trabajados) basta:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Regresión polinómica (transformación polinómica de las features)", fontsize=14, fontweight="bold")
### Training
axes1=plt.subplot(1,2,1, title="Training")
# Dibujamos los puntos (aquí el orden da igual):
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)
# Ordenamos los datos por el eje x con el método
# .sort_values() de los DataFrames:
dataset_poly2_training_sorted = dataset_poly2_training.sort_values(by="años trabajados")
# Dibujamos la curva de regresión polinómica:
plt.plot(dataset_poly2_training_sorted["años trabajados"],
lr_2.predict(dataset_poly2_training_sorted[["años trabajados","años trabajados ^2"]]),
c="red")
### Testing
axes2=plt.subplot(1,2,2, title="Testing")
# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
# Y también ordenamos los de testing:
dataset_poly2_testing_sorted = dataset_poly2_testing.sort_values(by="años trabajados")
# Dibujamos la curva de regresión polinómica:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
lr_2.predict(dataset_poly2_testing_sorted[["años trabajados","años trabajados ^2"]]),
c="red")
pass
No dibujamos el cuadrado de los años trabajados. Para temas de visualización, queda mejor si solo utilizamos las features originales (realmente, hemos usado un modelo lineal sobre un dataset transformado de forma no lineal).
Vemos que ajusta mucho mejor. Podríamos probar también con polinomios de grado 3, 4, 5... Pero si elevamos mucho el grado del polinomio, acabaríamos sobreajustando el modelo de nuevo.
Probemos otros regresores distintos.
6.D. Otros regresores¶
Sklearn incluye más regresores aparte de la regresión lineal, su posible expansión polinómica y sus regularizaciones. Veamos algunos:
6.D.a. Árboles regresores¶
sklearn.tree.DecisionTreeRegressor nos permite usar árboles para problemas de regresión. Los árboles son modelos no lineales, así que podemos olvidarnos de tema de transformaciones no lineales del dataset:
from sklearn.tree import DecisionTreeRegressor
# Formulamos árbol:
tree_reg = DecisionTreeRegressor(criterion="mse", # Al ser árboles regresores, no usan gini o entropía, sino MSE
max_depth=4,
min_samples_split=5,
min_samples_leaf=5
)
# Entrenamos:
tree_reg.fit(X=dataset_regresion_training[["años trabajados"]],
y=dataset_regresion_training["salario bruto"])
# Medimos:
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
y_pred=tree_reg.predict(dataset_regresion_testing[["años trabajados"]]))
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " + str(rmse))
Dibujamos:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Árbol regresor (profundidad máxima de 4)", fontsize=14, fontweight="bold")
# De nuevo, al no ser una línea recta
# la solución, debemos utilizar el dataset
# ordenado en el eje x de menor a mayor
# (ya lo tenemos ordenado de antes):
### Training
axes1=plt.subplot(1,2,1, title="Training")
# Dibujamos los puntos
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)
# Dibujamos la regresión de árbol (con el dataset ordenador por la x):
plt.plot(dataset_poly2_training_sorted["años trabajados"],
tree_reg.predict(dataset_poly2_training_sorted[["años trabajados"]]),
c="red")
### Testing
axes2=plt.subplot(1,2,2, title="Testing")
# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
# Dibujamos la regresión de árbol:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
tree_reg.predict(dataset_poly2_testing_sorted[["años trabajados"]]),
c="red")
pass
Las reglas de decisión y splits de este árbol se pueden dibujar de la misma forma que las de los árboles clasificadores (cualquier modelo que proviene de sklearn.tree puede ser dibujado con el mismo procedimiento).
6.D.b. Ensembles¶
Al igual que existen ensembles de clasificadores, existen ensembles de regresores. Para los ensembles de bagging la técnica suele ser hacer la media de las predicciones de cada modelo lanzado. En boosting, se intenta ir corrigiendo el $MSE$ de los modelos anteriores con los nuevos.
Probemos un Random Forest regresor sklearn.ensemble.RandomForestRegressor:
from sklearn.ensemble import RandomForestRegressor
rf_reg = RandomForestRegressor(n_estimators=20, # 20 árboles
criterion="mse",
max_depth=3,
min_samples_split=10,
min_samples_leaf=5,
bootstrap=True)
# Entrenamos:
rf_reg.fit(X=dataset_regresion_training[["años trabajados"]],
y=dataset_regresion_training["salario bruto"])
# Medimos:
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
y_pred=rf_reg.predict(dataset_regresion_testing[["años trabajados"]]))
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " + str(rmse))
Dibujamos:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Random Forests (regresor)", fontsize=14, fontweight="bold")
### Training
axes1=plt.subplot(1,2,1, title="Training")
# Dibujamos los puntos
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)
# Dibujamos la regresión:
plt.plot(dataset_poly2_training_sorted["años trabajados"],
rf_reg.predict(dataset_poly2_training_sorted[["años trabajados"]]),
c="red")
### Testing
axes2=plt.subplot(1,2,2, title="Testing")
# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
# Dibujamos la regresión:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
rf_reg.predict(dataset_poly2_testing_sorted[["años trabajados"]]),
c="red")
pass
6.D.c. Regresión por Nearest Neighbors¶
Su funcionamiento es sencillo: mira los vecinos cercanos a los puntos a predecir, y hace la media de los valores de dichos vecinos. Al igual que su homólogo clasificador, es un modelo no lineal.
Se encuentra en sklearn.neighbors.KNeighborsRegressor:
from sklearn.neighbors import KNeighborsRegressor
knn_reg = KNeighborsRegressor(n_neighbors=5, # 5-NN
metric="euclidean")
# Entrenamos:
knn_reg.fit(X=dataset_regresion_training[["años trabajados"]],
y=dataset_regresion_training["salario bruto"])
# Medimos:
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
y_pred=knn_reg.predict(dataset_regresion_testing[["años trabajados"]]))
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " + str(rmse))
figura = plt.figure(figsize=(16,7))
figura.suptitle("5-NN (regresor)", fontsize=14, fontweight="bold")
### Training
axes1=plt.subplot(1,2,1, title="Training")
# Dibujamos los puntos
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)
# Dibujamos la regresión:
plt.plot(dataset_poly2_training_sorted["años trabajados"],
knn_reg.predict(dataset_poly2_training_sorted[["años trabajados"]]),
c="red") # En rojo se ve mejor
### Testing
axes2=plt.subplot(1,2,2, title="Testing")
# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
# Dibujamos la regresión:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
knn_reg.predict(dataset_poly2_testing_sorted[["años trabajados"]]),
c="red") # En rojo se ve mejor
pass
6.D.d. Support Vector Regression¶
Aunque no tan famosos como sus homólogos clasificadores, los SVRs (Support Vector Regressors) también funcionan y generalizan muy bien. El que vamos a usar está en sklearn.svm.SVR:
from sklearn.svm import SVR
svr = SVR(C=0.5,
kernel="poly", # kernel polinómico
degree=2, # de segundo grado
)
# Entrenamos:
svr.fit(X=dataset_regresion_training[["años trabajados"]],
y=dataset_regresion_training["salario bruto"])
# Medimos:
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
y_pred=svr.predict(dataset_regresion_testing[["años trabajados"]]))
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " + str(rmse))
figura = plt.figure(figsize=(16,7))
figura.suptitle("SVR (kernel polinómico de grado 2)", fontsize=14, fontweight="bold")
### Training
axes1=plt.subplot(1,2,1, title="Training")
# Dibujamos los puntos
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)
# Dibujamos la regresión:
plt.plot(dataset_poly2_training_sorted["años trabajados"],
svr.predict(dataset_poly2_training_sorted[["años trabajados"]]),
c="red") # En rojo se ve mejor
### Testing
axes2=plt.subplot(1,2,2, title="Testing")
# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
# Dibujamos la regresión:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
svr.predict(dataset_poly2_testing_sorted[["años trabajados"]]),
c="red") # En rojo se ve mejor
pass
Hay muchos más. Cada vez que veas un modelo con la coletilla Regressor... Ya sabes :)
7. Clustering¶
A diferencia de la clasificación y la regresión, los modelos de clustering entran dentro de lo que se llama aprendizaje no supervisado. Tanto en clasificación como en regresión, entrenamos los modelos con datos de los que sabemos el valor real a predecir:
- En el caso de clasificación los hemos llamado labels o clases (calvo/no-calvo en nuestro ejemplo).
- En el caso de regresión los hemos llamado targets (el salario bruto anual).
En el aprendizaje no supervisado no es así: no tenemos un valor real a predecir: bien porque no está disponible en los datos, o bien porque no existe como tal. Por eso en aprendizaje no supervisado no se suele hablar de modelos precisos o imprecisos, sino de modelos útiles o apropiados para lo que queremos conseguir.
Un ejemplo de aprendizaje no supervisado son los modelos de clustering. Su principio es sencillo: buscar entre nuestros datos, y agruparlos de forma que cada grupo sea lo más homogéneo posible, y lo más diferente posible del resto de grupos.
A mí me gusta dividir los algoritmos/modelos de clustering en dos grandes grupos:
- Clusterings basados en distancia: al fin y al cabo, nuestras observaciones se pueden representar como puntos en un espacio dimensional, donde las dimensiones son tantas como el número de features que tenemos. Los modelos de clustering basados en distancia utilizan una métrica de distancia (como puede ser la euclídea, la Manhattan o la Mahalanobis) para separar los distintos grupos o clusters.
- Clusterings basados en densidad: computan los grupos no teniendo en cuenta solo la distancia absoluta entre los puntos, sino detectando "zonas" donde la densidad de puntos/observaciones es mayor.
Sklearn tiene modelos de clustering de ambos tipos. Preparemos un par de datasets, y pongámonos manos a la obra:
7.A. Los datos¶
Utilizaremos dos datasets: por un lado uno sintético creado con sklearn.datasets.make_blobs(); y por otro lado otro preparado para probar algoritmos de clustering. Comencemos por el de make_blobs:
dataset = make_blobs(n_samples=1350,
n_features=2,
centers=6,
cluster_std=1.5,
random_state=11)
Recuerdo que make_blobs nos devuelve una tupla con dos arrays de Numpy: uno con las features, y otro con las labels. Vamos a ignorar completamente las labels, al estar en aprendizaje no supervisado. Vamos a cargar simplemente las features en un DataFrame:
dataset_clust = pd.DataFrame(dataset[0])
# Nombramos columnas. Paso de inventarme historias:
# les llamo x e y a las features:
dataset_clust.columns = ["x","y"]
dataset_clust[:5]
Dibujamos:
dataset_clust.plot(kind="scatter",
x="x",
y="y",
figsize=(8,6))
pass
El random_state o seed del generador de números aleatorios ha sido elegido a propósito, para que los dos blobs del centro parezcan más o menos uno solo. Pero no me gusta: demasiado sencillo.
Vamos a generar algo de ruido, metiendo un blob grande y con gran desviación típica, y lo añadimos al dataset, a modo de "datos huérfanos", sin cluster ovbio:
ruido = make_blobs(n_samples=550,
n_features=2,
centers=1,
cluster_std=8,
center_box=[0,0],
random_state=19)
# Lo metemos en un DataFrame:
ruido_df = pd.DataFrame(ruido[0])
ruido_df.columns = ["x","y"]
# Y lo concatenamos con dataset_clust:
dataset_clust_1 = pd.concat([dataset_clust, ruido_df])
Y dibujamos:
dataset_clust_1.plot(kind="scatter",
x="x",
y="y",
figsize=(8,6))
pass
Mucho mejor. Precioso.
Por otro lado, cargamos el dataset preparado (que resulta estar en formato csv). Estos datos son también sintéticos, cedidos amablemente por los creadores de la implementación en Python del HDBSCAN:
dataset_clust_2 = pd.read_csv("clustering_data.csv", sep=";")
dataset_clust_2[:5]
Y lo dibujamos:
dataset_clust_2.plot(kind="scatter",
x="x",
y="y",
figsize=(8,6))
pass
Mola porque tiene clusters no esféricos, así como ruido.
Ya podemos empezar a lanzar modelos.
7.B. Algoritmos de clustering¶
Aquí muestro algunos de los que incluye sklearn (como digo siempre, hay más de los listados aquí):
| Algoritmo | Módulo |
|---|---|
| K-Means | cluster.KMeans |
| Affinity propagation | cluster.AffinityPropagation |
| Mean-shift | cluster.MeanShif |
| Spectral clustering | cluster.SpectralClustering |
| Agglomerative clustering (jerárquico) | cluster.AgglomerativeClustering |
| DBSCAN | cluster.DBSCAN |
| Gaussian mixtures | mixture.GMM |
| Birch | cluster.Birch |
Probemos algunos de estos con nuestros dos datasets.
7.B.a. K-Means¶
El K-Means es el algoritmo de clustering más conocido. Sencillo de entender, basado en distancia euclídea y computacionalmente no muy intensivo (si bien se trata de un problema que se resuelve a través de métodos heurísticos, el Algoritmo de Lloyd hace un buen trabajo convergiendo para el K-Means en la mayoría de los casos). A lo largo de estas iteraciones, los centroides de los clusters se van moviendo, siendo la media de los puntos de cada cluster en la iteración anterior. Y de ahí el nombre K-means (o K-medias).
Usémoslo en nuestros dos datasets con sklearn.cluster.KMeans:
from sklearn.cluster import KMeans
km1 = KMeans(n_clusters=6, # Pega del KMeans: nos obliga a elegir número de clusters.
max_iter=400)
# Entrenamos con el dataset_clust_1:
km1.fit(dataset_clust_1) # OJO: cada vez que se entrena, los clusters "no coinciden",
# es decir: lo que en un entrenamiento es el cluster número 3,
# en otro entrenamiento el cluster será el mismo (que contiene
# más o menos los mismos puntos), pero igual pasa a ser el número
# 2, o 4, o el que sea...
# Ocurre en todos los modelos de clustering no-determinísticos.
# Y predecimos a qué cluster pertenece cada punto
# que se representan con números, comenzando por el 0
# hasta n_clusters - 1:
result_kmeans_dataset1 = km1.predict(dataset_clust_1)
result_kmeans_dataset1
Podemos utilizar estas etiquetas para pintar a qué cluster pertenece cada punto:
dataset_clust_1.plot(kind="scatter",
x="x",
y="y",
c=km1.predict(dataset_clust_1),
cmap=color_maps.Accent,
figsize=(9,7))
pass
Y para el otro dataset (donde asumimos que también hay 6 clusters):
km2 = KMeans(n_clusters=6,
max_iter=600)
# Entrenamos con el dataset_clust_2:
km2.fit(dataset_clust_2)
# Y dibujamos:
dataset_clust_2.plot(kind="scatter",
x="x",
y="y",
c=km2.predict(dataset_clust_2),
cmap=color_maps.Accent,
figsize=(9,7))
pass
Rápidamente podemos ver los problemas del K-Means:
- Necesitamos decirle el número de clusters a crear. Existen técnicas para intentar inferir el número óptimo de clusters, pero en la práctica ninguna es maravillosa. Cierto es que en ocasiones el número de clusters es una imposición (por ejemplo, en una campaña de segmentación de marketing donde hay presupuesto para tratar de forma diferenciada a X grupos de clientes), pero en el resto de casos siempre es complicado...
- Los clusters generados siempre son esféricos (o circulares en el caso de dos dimensiones).
- Todos los puntos pertenecerán a algún cluster: el modelo no descarta ningún punto; por muy alejado que esté del resto. En casos muy extremos, puede hacer que un punto (o unos cuantos) muy alejados formen su propio cluster alejado; pero ningún punto es descartado.
7.B.b. Gaussian Mixture Model¶
El Gaussian Mixture Model (abreviado como GMM) es un modelo probabilístico que es usado comúnmente como modelo de clustering. De hecho, se trata de la generalización del K-means, o dicho de otro modo: el K-Means es un caso muy específico de GMM.
Probemos el sklearn.mixture.GaussianMixture (sí: no está en sklearn.cluster, sino en sklearn.mixture):
from sklearn.mixture import GaussianMixture
gmm1 = GaussianMixture(n_components=6, # En GMMs se habla de componentes en vez de clusters
covariance_type="diag", # Tipo de parámetros de covarianza a usar
max_iter=250) # Número máximo de iteraciones
# Entrenamos con el dataset_clust_1:
gmm1.fit(dataset_clust_1)
# Y dibujamos:
dataset_clust_1.plot(kind="scatter",
x="x",
y="y",
c=gmm1.predict(dataset_clust_1),
cmap=color_maps.Accent_r,
figsize=(9,7))
pass
GMM ha determinado el ruido como un cluster... A ver con el otro dataset:
gmm2 = GaussianMixture(n_components=6,
covariance_type="diag",
max_iter=500)
# Entrenamos con el dataset_clust_2:
gmm2.fit(dataset_clust_2)
# Y dibujamos:
dataset_clust_2.plot(kind="scatter",
x="x",
y="y",
c=gmm2.predict(dataset_clust_2),
cmap=color_maps.Accent,
figsize=(9,7))
pass
Bastante similar al K-Means en este caso. Consideraciones del GMM:
- Necesitamos definir el número de clusters a crear, igual que en el K-Means.
- Los clustes no tienen por qué ser esféricos: el K-Means contempla únicamente la media de los puntos para determinar los clusters (el vector de medias), mientras que el GMM contempla también la matriz de varianzas y covarianzas de los datos. Es esto lo que le permite la creación de clusters no esféricos.
- Todos los puntos pertenencen a algún cluster; no se descartan puntos.
- A diferencia del K-Means, la salida real del modelo no es una etiqueta de a qué componente/cluster pertenece cada punto; sino una distribución de probabilidad, donde se muestra la probabilidad de que cada punto pertenezca a cada componente/cluster. Eso sí: automáticamente, sklearn (y todas las bibliotecas que implementan GMM) nos da el componente/cluster al que es más probable que pertenezca el punto. No obstante, el método
.predict_proba()nos ofrece esta distribución de probabilidad (¡otros modelos, incluso clasificadores, también lo ofrecen!)
7.B.c. Affinity Propagation¶
El Affinity Propagation es un modelo que se basa en grafos y que hace a cada punto votar a su cluster o exemplar favorito. A diferencia de en el K-Means y GMM, los centroides (o centros de los clusters) en el Affinity Propagation son puntos reales del propio dataset, los cuales son llamados exemplars. A partir de ahí, funciona de forma bastante similar al K-Means. Eso sí: la formulación del algoritmo permite utilizarlo con métricas no distanciales, esto es, que la métrica de similaridad/no-similaridad utilizada no tiene por qué satisfacer las cualidades de las distancias.
Utilicemos pues sklearn.cluster.AffinityPropagation:
from sklearn.cluster import AffinityPropagation
ap1 = AffinityPropagation(damping=0.9, # Damping factor, >=0.5 y <1. A más grande, menos exemplars.
affinity="euclidean", # Métrica de distancia euclídea
preference=-1200.0, # Preferencias para cada punto (si es array),
# o generales(si es float). Por lo general,
# recomiendo un número negativo, más grande o
# más pequeño dependiendo de los datos. Si no,
# nos sacará muchísimos clusters...
max_iter=3000) # Máximo de iteraciones
# Entrenamos con el dataset_clust_1:
ap1.fit(dataset_clust_1)
# Y dibujamos:
dataset_clust_1.plot(kind="scatter",
x="x",
y="y",
c=ap1.predict(dataset_clust_1),
cmap=color_maps.Accent,
figsize=(9,7))
pass
Sobre el otro dataset:
ap2 = AffinityPropagation(damping=0.95,
affinity="euclidean",
preference=-3.0,
max_iter=3000)
# Entrenamos con el dataset_clust_2:
ap2.fit(dataset_clust_2)
# Y dibujamos:
dataset_clust_2.plot(kind="scatter",
x="x",
y="y",
c=ap2.predict(dataset_clust_2),
cmap=color_maps.Accent,
figsize=(9,7))
pass
Sobre el Affinity Propagation:
- No necesitamos decirle el número de clusters o exemplars: el propio modelo los determina a partir de la medida de preferencia y el dampening factor.
- Asume clusters esféricos.
- Todos los puntos pertenecen a algún cluster; no se descartan puntos.
- Es lento: tiene que determinar qué puntos van a ser exemplars, tanto en cantidad como en localización, y luego ya determinar a qué exemplar pertenece cada uno del resto de puntos.
7.B.d. Mean Shift¶
El Mean Shift es un modelo de clustering bastante potente. Básicamente el modelo asume que existe cierta función de densidad (probabilística) que explica la distribución de los datos, e intenta colocar centroides para maximizar la verosimilitud en función de los datos. Realiza dicha aproximación utilizando técnicas de KDE (Kernel Density Estimation).
El modelo se encuentra en sklearn.cluster.MeanShift:
from sklearn.cluster import MeanShift
ms1 = MeanShift(bandwidth=5.5, # El "ancho de banda" del kernel utilizado para la estimación
cluster_all=False) # Podemos elegir que "clusterice" todos los puntos, o no
## Para los modelos de clustering sobre todo:
# existe un método llamado f.it_predict(), que
# entrena y predice sobre los datos dados automáticamente.
# No es ideal cuando tenemos training y testing, pero sí
# en modelos no supervisados. En vez de hacer .fit() y
# luego .predict() sobre los mismos datos, hacemos
# .fit_predict() en un solo método:
dataset_clust_1.plot(kind="scatter",
x="x",
y="y",
c=ms1.fit_predict(dataset_clust_1),
cmap=color_maps.Accent_r, # Con el _r revertimos el orden de colores,
figsize=(9,7)) # lo cual en este caso nos pinta los -1 de color gris oscuro.
pass
¿Ves los puntos de color gris oscuro? Son los puntos descartados, que no pertenecen a ningún cluster (son etiquietados con el número $-1$ a modo de label en el resultado). Veámoslo en el otro dataset:
ms2 = MeanShift(bandwidth=0.17,
cluster_all=False)
# Entrenamos y dibujamos:
dataset_clust_2.plot(kind="scatter",
x="x",
y="y",
c=ms2.fit_predict(X=dataset_clust_2),
cmap=color_maps.Accent_r,
figsize=(9,7))
pass
Misma historia: los puntos gris oscuro han sido descartados.
Con el Mean Shift:
- No tenemos que determinar el número de clusters. Pero al igual que con el Affinity Propagation, tenemos que "tener cuidado" con los parámetros (el
bandwidthen este caso, principalmente) para obtener una solución razonable. - Los clusters resultantes son esféricos.
- Es capaz de descartar puntos, que considera que no pertenecen a ningún cluster (de forma opcional).
7.B.e. DBSCAN¶
El DBSCAN es un algoritmo de clustering basado en densidad: asume que hay clusters en regiones donde la densidad de los puntos es elevada. DBSCAN genera una transformación del espacio dimensional basándose en la densidad de los datos: los puntos en zonas densas son separados, mientras que los de zonas menos densas son apartados de forma iterativa.
Utilicemos sklearn.cluster.DBSCAN:
from sklearn.cluster import DBSCAN
dbs1 = DBSCAN(eps=0.95, # Parámetro epsilon, que cambia la "consideración de densidad" del algoritmo
min_samples=15) # Número de puntos vecinos necesarios para que un punto sea "core" de densidad
# Entrenamos y dibujamos:
dataset_clust_1.plot(kind="scatter",
x="x",
y="y",
c=dbs1.fit_predict(dataset_clust_1),
cmap=color_maps.Accent_r,
figsize=(9,7))
pass
De nuevo: puntos grises, descartados. Probemos con el dataset complicado:
dbs2 = DBSCAN(eps=0.035,
min_samples=20)
# Entrenamos y dibujamos:
dataset_clust_2.plot(kind="scatter",
x="x",
y="y",
c=dbs2.fit_predict(dataset_clust_2),
cmap=color_maps.Accent_r,
figsize=(9,7))
pass
Sobre el DBSCAN:
- "Tocando" los parámetros epsilon y el número mínimo de puntos por core de densidad (no deberíamos hablar de centroides en DBSCAN), el modelo inferirá el número de clusters necesarios.
- Los clusters resultantes no tienen que ser esféricos (a la vista está).
- Descarta zonas de puntos donde la densidad es baja (descarta puntos por tanto).
Recuerdo: en el mundo del aprendizaje no supervisado (dentro del cual se incluye el clustering) NO hay una solución correcta; sino soluciones más convenientes que otras para un determinado propósito.
Recomiendo también ver (aparte del resto de modelos de clustering de sklearn) HDBSCAN, que es un modelo muy similar al DBSCAN. La diferencia es que HDBSCAN permite tener clusters densos, pero de densidades diferentes entre ellos.
No está incluído en sklearn, sino como biblioteca aparte. Sin embargo, se utiliza exactamente igual que los modelos de clustering de sklearn. E instalarlo es tan sencillo como:
$ pip install hdbscan
Nota a futuro para ti para siempre: sé que es tentador. Sé que te lo van a pedir en el trabajo. Sé que es... incluso divertido. Pero por favor: NO INCLUYAS NUNCA VARIABLES CATEGÓRICAS/DICOTÓMICAS EN UN MODELO DE CLUSTERING. La razón es sencilla; pongámonos en el caso de variables dicotómicas (que son o $1$ o $0$): en un dataset compuesto por dos features dicotómicas, de forma que cada observación puede ser $0|0$, $0|1$, $1|0$ o $1|1$... ¿Cuántos clusters hay? Pues 4, evidentemente: las 4 combinaciones. Introducir variables no numéricas en un modelo de clustering solo consigue una cosa: contaminar el modelo. Si el objetivo es hacer clustering y tienes variables numéricas y categóricas/dicotómicas, haz primero clustering sobre las numéricas, y luego computa todas las posibles combinaciones de variables categóricas/dicotómicas para cada cluster generado. Esto funciona, y es lo razonable.
Y con esto concluimos la parte de clustering. Seguimos con el feature engineering.
8. Feature Engineering¶
Lanzar modelos y obtener resultados es la parte sencilla en machine learning. El verdadero reto muchas veces es obtener features significativas y útiles para nuestros modelos.
Cuando te enfrentas a un problema de data science de primeras, normalmente comienzas por obtener los datos "en bruto": datos transaccionales mes a mes, comportamientos de cliente, imágenes de distintos tamaños y con diferentes perspectivas, etcétera.
Muchas veces el primer paso es generar variables derivadas a partir de estos datos en bruto (diferencias mes a mes, combinación de variables de crédito/débito, valores absolutos de balances, normalización del perfil RGB de fotografías...). Este tipo de transformaciones de los datos están muy vinculadas al sector/negocio en el que nos movemos, más que al punto de vista matemático.
Luego, es frecuente "filtrar" features: aquellas que tienen gran porcentaje de nulos, elevada correlación con otras, o valores constantes (que no permiten discriminar de ninguna manera).
A partir de ahí, hay tres formas principales de preparar los datasets para nuestros modelos:
- Incluir todas las variables (tanto en crudo o derivadas) válidas, y utilizar alguno de los modelos que "solo prestan atención" a las variables significativas. Ejemplos de estos modelos son aquellos que utilizan regularización $L1$ (como la regresión LASSO), los árboles (o ensembles de los mismos) y las redes.
- Utilizar modelos de reducción de dimensionalidad: si tenemos cientos (o miles, o incluso millones) de features, la reducción de dimensionalidad permite generar combinaciones de todas ellas, generando un dataset de variables virtuales mucho más compacto, más compatible con la mayoría de clasificadores y regresores. Los modelos de reducción de dimensionalidad son modelos de aprendizaje no supervisado (igual que los de clustering).
- Seleccionar variables, utilizando los propios resultados arrojados por los modelos. De esta forma podemos ir decidiendo si añadir o quitar más o menos variables, en función de si los resultados van mejorando o empeorando.
Realmente no sueles limitarte a tomar uno de estos tres caminos (¡sobre todo si es el primero!). Lo normal es usar una combinación de estas tres técnicas.
Nota: los modelos de reducción de dimensionalidad reducen drásticamente la explicabilidad del modelo. A veces esto no es relevante; pero otras veces sí lo es. Depende mucho del negocio y del valor que se quiere que aporte el machine learning. Si que un modelo sea opaco o difícil de entender no es problema, entonces bien. Pero en otros ámbitos (por ejemplo el bancario, donde si no ofreces un crédito a una persona debes decirle por qué) la explicabilidad y transparencia del modelo es esencial. En estos casos, la reducción de dimensionalidad no suele ser una opción viable, y el camino a elegir suele ser la selección de variables.
En este apartado vamos a explorar la reducción de dimensionalidad y la selección de variables.
8.A. Reducción de dimensionalidad¶
Si bien no hay una definición formal que englobe a los algoritmos de reducción de dimensionaliad (al menos que yo conozca), dichos algoritmos son fácilmente reconocibles: a partir de una serie de variables o features, genera otras completamente diferentes e independientes, intentanto extraer la relevancia y cualidades determinantes de las originales. El resultado es una menor cantidad de variables no originales (conocidas a veces como variables virtuales), pero muy significativas.
Estas variables virtuales generadas no tienen significado por sí mismas; son solo features generadas y que se asume que son significativas y útiles.
El modelo más famoso y utilizado en reducción de dimensionalidad es el PCA (Principal Component Analysis) y sus variantes, disponibles en sklearn.decomposition.
8.A.a. Principal Component Analysis¶
PCA o Principal Component Analysis transforma las variables o features que pasemos a dicho modelo en unas cuantas features virtuales (menos que las originales). Podemos pensar en PCA como en un pliegue de las variables en la dirección de la máxima varianza posible. De esta forma, las features virtuales generadas tendrán la mayor varianza posible; lo cual suele dar mayor poder discriminador al análisis que hagamos a posteriori. Supongamos que tenemos el siguiente dataset, con dos features $x_1$ y $x_2$:
Y queremos obtener una única feature virtual a partir de esas dos. PCA va a intentar plegar los puntos en la dirección en la cual la varianza de la nueva feature resultante es la máxima posible.
Supongamos que realizamos el siguiente pliegue:
La línea azul muestra el "nuevo eje de coordenadas" de la nueva variable generada. Como solo queríamos una, solo habrá un eje. No obstante, podemos ver que la distribución de esta nueva variable generada (los puntos azules) no tiene demasiada varianza; todos los puntos están muy cerca unos de otros. Esto es porque no hemos plegado en la dirección de la máxima varianza de las features originales.
Si de verdad utilizamos PCA, el resultado sería algo así:
Ahora sí hemos plegado las dos features originales en la dirección de la máxima varianza posible. Los puntos azules tienen ahora mucha más varianza. Intuitivamente podemos ver que esta nueva feature virtual generada captura la mayoría de la varianza de las features originales.
Y esto es exactamente lo que hace PCA: intentar generar variables virtuales que aporten toda la información posible sobre las variables originales, pero en una representación más compacta (en este caso pasando de tener dos variables a solo una). Esto se llama utilizar la top principal component del dataset original.
En un caso donde tuviéramos originalmente 300 variables y quisiéramos quedarnos con 5, usaríamos PCA para quedarnos con las top 5 principal components del dataset.
Vamos a utilizar sklearn.decomposition.PCA sobre el dataset de clasificación que habíamos generado en el apartado de regresión. Era este (ahora mismo las labels nos dan igual):
dataset_clasificacion.plot(kind="scatter",
x="hormona_a",
y="hormona_b",
edgecolors="black",
linewidths=1,
s=30,
figsize=(8,6))
pass
Apliquemos PCA para obtener la top (1) principal component del dataset:
from sklearn.decomposition import PCA
pca = PCA(n_components=1) # La top 1 principal component
# PCA tiene los métodos .fit() para entrenar,
# .transform() para transformar el dataset y obtener
# nuestras top components, y .fit_transform() para
# hacer ambas en un solo paso:
top_principal_component = pca.fit_transform(dataset_clasificacion[["hormona_a", "hormona_b"]])
top_principal_component es un array con la nueva variable virtual generada (la top principal component). Vamos a meterla en un DataFrame, junto con las dos features originales:
dataset_clasificacion_con_pca = dataset_clasificacion.copy() # Así podemos copiar un DataFrame
dataset_clasificacion_con_pca["top 1 principal component"] = top_principal_component
# Reordenamos columnas para que quede más mono:
dataset_clasificacion_con_pca = dataset_clasificacion_con_pca[["hormona_a", "hormona_b", "top 1 principal component", "label"]]
dataset_clasificacion_con_pca[:10]
Vamos a ver la varianza:
dataset_clasificacion_con_pca.var() # Método .var() del DataFrame
Podemos ver que la top principal component intenta capturar toda la varianza posible de hormona_a y hormona_b, dando como resultado una variable con más varianza que las dos originales. Ahora podríamos utilizar esta principal component como variable, en vez de las dos originales; de esta forma, tendríamos un problema con una dimensionalidad menor (1-D en vez de 2-D).
PCA también se utiliza mucho para visualización: imagínate que tenemos un dataset con 500 features, y queremos dibujar un scatterplot... Pues es imposible: solo podemos representar hasta 3 dimensiones de forma fácilmente entendible. Pues bien, podríamos utilizar PCA para obtener las top 3 principal components del dataset, y dibujar el scatterplot de los datos en función de estas 3 dimensiones creadas.
8.B. Selección de variables¶
Los modelos de reducción de dimensionalidad como PCA no deben ser utilizados a la ligera. Como hemos comentado antes, las variables generadas son virtuales, y carecen de significado como tal. En determinados ámbitos esto es aceptable (por ejemplo en neurociencia, donde las variables originales son zonas del cerebro; y hay miles de estas). Pero en otros ámbitos más estándar (marketing, banca y seguros, sector media...) normalmente queremos modelar con variables reales.
En estos casos, si queremos trabajar con menos variables, suele tocar seleccionar.
Para este ejemplo vamos a utilizar el dataset de los precios de las casas de Boston:
boston_training_set[:5]
Donde tenemos 13 features más nuestro target.
sklearn.feature_selection incluye varias técnicas para hacer selección de variables. Vamos a probar algunas.
Pero primero vamos a usar el poder de los DataFrames jerárquicos para hacer nuestras vida más sencilla. Vamos a generar una jerarquía de columnas que nos permita separar las features del target fácilmente:
columnas_jerarquicas = pd.MultiIndex.from_tuples([("features", "CRIM"),
("features", "ZN"),
("features","INDUS"),
("features","CHAS"),
("features","NOX"),
("features","RM"),
("features", "AGE"),
("features", "DIS"),
("features", "RAD"),
("features", "TAX"),
("features", "PTRATIO"),
("features", "B"),
("features", "LSTAT"),
("target", "Precio")])
# Y las aplicamos:
boston_training_set.columns = columnas_jerarquicas
boston_training_set[:5]
Mejor. Vamos a ello.
8.B.a. Filtro de variables por varianza¶
sklearn.feature_selection.VarianceThreshold nos elimina automáticamente las variables con baja varianza. La intuición es sencilla: variables que no varían mucho a lo largo del dataset, son poco informativas:
from sklearn.feature_selection import VarianceThreshold
filtro_var = VarianceThreshold(threshold=5.0) # Todas las features con menor varianza que
# este número, serán filtradas.
# Lo "entrenamos" (no es un entrenamiento como tal,
# pero configura el selector de variables para nuestro
# dataset):
filtro_var.fit(boston_training_set["features"])
# Y transformamos nuestro dataset, efectivamente descartando
# variables:
boston_var_filtered = filtro_var.transform(boston_training_set["features"])
# La primera fila tras la transformación:
boston_var_filtered[0]
Ha dejado nuestras 13 features en 8... Pero resulta complicado (o tedioso) saber cuáles ha filtrado, puesto que hemos perdido los nombres de las columnas.
Esta función permite recuperarlas:
def recuperar_columnas(dataframe_original, input_array):
resultado = pd.DataFrame()
input_dataframe = pd.DataFrame(input_array)
input_dataframe.columns = [str(col) for col in input_dataframe.columns]
for col_orig in dataframe_original.columns:
for col_reduc in input_dataframe.columns:
if np.all(dataframe_original[col_orig].values==input_dataframe[col_reduc].values):
resultado[col_orig] = dataframe_original[col_orig]
else:
continue
# Si es un DataFrame jerárquico:
if type(resultado.columns[0])==tuple:
multiIndex = pd.MultiIndex.from_tuples(resultado.columns)
resultado.columns = multiIndex
return resultado
La aplicamos:
boston_var_filtered = recuperar_columnas(boston_training_set, boston_var_filtered)
# Y volvemos a añadir el target:
boston_var_filtered["target","Precio"] = boston_training_set["target","Precio"]
boston_var_filtered[:5]
El método de filtrar variables por varianza está bien, pero se basa bastante el el criterio personal de: ¿qué es una varianza alta?
8.B.b. Selección de variables univariante¶
sklearn presenta varias funciones para seleccionar variables en función de pruebas realizadas a las variables una por una (es decir, sin medir interacciones de ningún modo entre ellas).
Uno de ellos es sklearn.feature_selection.SelectKBest, que selecciona $K$ variables (nosotros elegimos $K$) a partir de una función de scoring que devuelva p-valores:
- Para regresión:
sklearn.feature_selection.f_regression - Para clasificación:
sklearn.feature_selection.chi2osklearn.feature_selection.f_classif
from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import f_regression
# Generamos la instancia de SelectKBest
kb = SelectKBest(score_func = f_regression, # Pasamos únicamente el método, sin argumentos
k=4 # Queremos quedarnos con las 4 variables mejores
)
Y ahora utilizamos .fit_transform por ejemplo, para realizar la selección:
boston_KBest = kb.fit_transform(X=boston_training_set["features"],
y=boston_training_set["target","Precio"])
# Y lo pasamos a DataFrame:
boston_KBest = recuperar_columnas(boston_training_set, boston_KBest)
# Y re-añadimos la columna de target:
boston_KBest["target","Precio"] = boston_training_set["target","Precio"]
boston_KBest[:5]
8.B.c. Uso de árboles u otros clasificadores y regresores para seleccionar variables¶
Como ya hemos visto, los árboles de decisión van realizando splits en función de las distintas variables para clasificar. Dichos splits comienzan por la variable más discriminante, y luego la siguiente, y luego la siguiente... Hasta que el árbol clasifica correctamente. Teniendo el orden de las variables por el que el árbol va haciendo los splits, tenemos el orden de importancia de las variables (desde el punto de vista de la opinión del árbol, que evidentemente no tiene por qué ser perfecto).
Los árboles no son los únicos modelos capaces de hacer esto. De hecho, cualquier modelo que tras ser entrenado con .fit() tenga los atributos .coef_ o .feature_importances_ puede ser utilizado para seleccionar variables.
La selección de variables acorde a un clasificador/regresor se puede hacer con sklearn.feature_selection.SelectFromModel:
from sklearn.feature_selection import SelectFromModel
# Creamos un arbol, que como el problema es de regresión...
arbol = DecisionTreeRegressor(criterion="mse", # Para que la métrica de criterio sea MSE
max_depth=6,
min_samples_split=5,
min_samples_leaf=5,
random_state=4
)
# Y se lo pasamos a una instancia de SelectFromModel:
sfm = SelectFromModel(estimator=arbol)
# Y usamos .fit_transform() para seleccionar:
boston_selectedFromModel = sfm.fit_transform(X=boston_training_set["features"],
y=boston_training_set["target","Precio"])
# Lo pasamos a DataFrame y añadimos columna de target:
boston_selectedFromModel = recuperar_columnas(boston_training_set, boston_selectedFromModel)
# Y re-añadimos la columna de target:
boston_selectedFromModel["target","Precio"] = boston_training_set["target","Precio"]
boston_selectedFromModel[:5]
Parece que el árbol ha decidido quedarse únicamente esas variables...
Ahora podríamos utilizar cualquier algoritmo (¡no necesariamente otro árbol!) para hacer nuestro modelo.
Recuerdo: el árbol construido y pasado a SelectFromModel solo se ha utilizado para seleccionar variables según su criterio.
8.B.d. Recursive Feature Elimination¶
La Recursive Feature Elimination es una de las técnicas de backward elimination para selección de variables. Comienzas con todas ellas, y pasas un clasificador (o regresor). Aquella variable que el modelo considera menos significativa es eliminada. Y así continuamente hasta que:
- Nos quedamos con el número de variables deseado.
- El modelo "decide" que no debe quitar ninguna más.
Para obtener el primer resultado (quedarnos con el número de features que queremos) utilizamos sklearn.feature_selection.RFE. Para obtener el segundo resultado se utiliza sklearn.feature_selection.RFECV: un Recursive Feature Elimination que además realiza validación cruzada para decidir cuántas variables debe dejar (más adelante explicaremos lo que es la validación cruzada).
Probemos los dos. Primero RFE:
from sklearn.feature_selection import RFE, RFECV
# Como modelo para realizar las RFEs vamos a usar un
# Support Vector Regressor:
mi_svr = SVR(kernel="linear")
# Y ahora instanciamos una RFE:
rfe = RFE(estimator = mi_svr,
n_features_to_select=7, # Quiero que se quede con 7
step=2 # Que vaya quitando variables de 2 en 2
)
# Y utilizamos .fit_transform() para ejecutar el RFE y que
# seleccione:
boston_selected_rfe = rfe.fit_transform(X=boston_training_set["features"],
y=boston_training_set["target","Precio"])
# Y a DataFrame(y metemos de nuevo el target):
boston_selected_rfe = recuperar_columnas(boston_training_set, boston_selected_rfe)
boston_selected_rfe["target","Precio"] = boston_training_set["target","Precio"]
boston_selected_rfe[:5]
Pidiéndole que nos dejara 7, el SVR ha decidido que deben ser esas 7. Pero... ¿y si no le decimos cuántas queremos, sino que le dejemos que él mismo decida? Para eso, utilizamos RFECV:
# Seguimos con el mismo SVR de antes, así que nada:
# creamos ya la instancia de RFECV:
rfecv = RFECV(estimator = mi_svr,
step=1, # Que vaya quitando de una en una
cv=10, # Número de "folds" de validación cruzada (lo explicamos luego)
scoring="neg_mean_squared_error" # Métrica de rendimiento. Un string de las de esta lista:
# http://goo.gl/bPGYAo
# o bien una función de las disponibles en sklearn.metrics
# (de las métricas de regresión en este caso, claro).
)
# Y utilizamos .fit_transform() para ejecutar el RFECV y que
# seleccione:
boston_selected_rfecv = rfecv.fit_transform(X=boston_training_set["features"],
y=boston_training_set["target","Precio"])
# Y a DataFrame(y metemos de nuevo el target):
boston_selected_rfecv = recuperar_columnas(boston_training_set, boston_selected_rfecv)
boston_selected_rfecv["target","Precio"] = boston_training_set["target","Precio"]
boston_selected_rfecv[:5]
Y decidió quedarse con esas.
Recomiendo prestar siempre muhcísima atención a la selección de variables. RFE y RFECV suelen funcionar muy bien. También es costumbre probarlos con diferentes clasificadores/regresores. Si todos ellos coinciden en determinada selección, es probablemente porque es la más adecuada.
Existen muchas técnicas más para la selección de variables. Una de ellas es la eliminación bidireccional o stepwise; que es similar a RFE, con la diferencia de que va eliminando pero también añadiendo variables, hasta probar todas las combinaciones posibles. Es computacionalmente muy cara (entrena el modelo para todas las posibles combinaciones de variables), pero resulta sencilla de implementar con un par de bucles for y los clasificadores y métricas de rendimiento que nos ofrece sklearn.
Dicho esto, pasemos a la validación y selección de modelos.
9. Validación y selección de modelos¶
Supongamos que ya hemos seleccionado variables para lanzar modelos. Conocemos cuáles hay, y sabemos más o menos cómo funcionan. Pero... ¿y todos los argumentos e hiperparámetros que toma cada modelo? ¿Cómo elijo eso?
La respuesta, tanto para seleccionar un modelo como sus hiperparámetros y configuraciones, es la validación.
La validación consiste en la aplicación de métricas para saber qué modelo (o variante de modelo) funciona mejor sobre datos que no ha visto antes. Parece que hace lo mismo que la separación en training y testing, ¿no?
Sí y no. Existen varios enfoques para la validación y testing de modelos. El más utilizado es el siguiente. Se parece a una especie de juego:
- Dividimos el dataset en dos conjuntos: el de training y testing. El conjunto de testing lo guardamos en una caja fuerte, y no vamos a volver a verlo ni utilizarlo hasta el final del proceso.
- Quedándonos con el conjunto de training, tenemos que hacer lo que podamos para elegir un modelo. Solo podemos elegir uno: tanto el modelo en sí como sus hiperparámetros y configuraciones.
- Cuando ya tengamos el modelo elegido (por ejemplo: regresión logística con regularización L2 y alpha de 0.1, con
fit_intercept=True), lo entrenamos con TODO nuestro training set. - Ahora ya, podemos abrir nuestra caja fuerte con el conjunto de testing. Los resultados de precisión obtenidos en este conjunto es lo que podemos esperar de nuestro modelo para el futuro.
Intuirás que lo complicado es el paso 2. Ahí es donde entra en juego la validación.
En validación, cogemos parte de nuestro training set y lo utilizamos para medir el rendimiento de nuestros modelos. Nos quedaremos con el que mejor rendimiento dé. De forma que lo que hemos ido haciendo en este Notebook ha sido una forma muy informal de validación; a pesar de haberla hecho con un conjunto de datos que hemos llamado testing (lo llamé testing para no meter al principio del todo el tema de la validación)
Repito: para ser léxicamente correctos: el conjunto de testing no debería verse ni utilizarse hasta que hayamos seleccionado nuestro modelo ganador mediante validación; y este conjunto de validación lo cogemos de una parte del conjunto de training.
9.A. Training, validación, testing¶
Hay distintas técnicas de validación. La más básica consiste en separar una parte del conjunto de training y llamarla conjunto de validación, de forma que:
- Entrenamos modelos con el conjunto de training.
- Comparamos, y seleccionamos el modelo que mejor rendimiento nos da en el conjunto de validación.
- Utilizamos training + validación para entrenar de nuevo el modelo ganador, sobre los dos conjuntos de datos a la vez.
- Medimos el rendimiento en el conjunto de testing; lo ponemos en una PPT y se lo entregamos a los inversores.
Para hacerlo, podemos utilizar sklearn.cross_validation.train_test_split. Lo hacemos una vez para separar en training y testing, y luego otra vez con el conjunto de training para separar en training de verdad y validación.
Lioso, lo sé. Este dibujito lo muestra de forma más sencilla:
Siempre surge la pregunta: ¿qué cantidad de datos dejar para cada conjunto? Existen varios compromisos:
- Si dejamos pocos datos para el test set y validation set, igual no son representativos del todo...
- Pero si cogemos muchos, los quitamos de training y nuestro modelo entrenará con menos datos...
Un criterio bastante utilizado es 60-20-20 (60% entrenamiento, 20% validación y 20% test). No funciona mal... Depende sobre todo de la cantidad de datos que tengamos. En general: cuantos menos datos, peor funciona este método (y todos; pero éste en especial).
Por suerte existe una forma alternativa, que permite que el conjunto de validación sea pequeño y grande a la vez... Y ahí entra en juego la validación cruzada.
9.B. Cross Validation¶
Cross Validation o validación cruzada son un conjunto de técnicas que hacen separaciones dinámicas y automáticas de conjuntos de training y validación.
Hay varias variantes de validación cruzada. Vamos a centrarnos en la más utilizada: K-Fold Cross Validation.
K-Fold Cross Validation consiste en que por cada modelo, no entrenamos una vez, sino un número $K$ de veces. Pongamos un valor de $K$ de 10. Pues bien, entrenaremos nuestro modelo 10 veces, cada vez dejando un 10% de los datos para validación; y cada vez siendo este 10% distinto. Una vez entrenado el modelo las 10 veces, sacamos el rendimiento del modelo haciendo la media de los resultados obtenidos ($MSE$, porcentaje bien clasificado, etc) en cada una de estas 10 veces. De esta forma, a nivel efectivo hemos entrenado con el 100% de los datos de training; y también hemos validado con el 100% de los datos de training (de nuevo: no tocamos el conjunto de testing).
Esta imagen (por cortesía de Neerav Basant) lo ilustra muy bien:
Podemos elegir el valor de $K$ que queramos. Lo normal es un valor entre 3 y 10. Cuanto mayor sea la $K$, más modelos serán entrenados, y más lento será el proceso (pero más robusto en teoría).
Exploremos ahora las distintas formas de hacer validación cruzada con sklearn:
9.B.a. Validación cruzada empotrada en la especificación de los modelos¶
Ciertos modelos de sklearn (principalmente los de sklearn.linear_model) tienen una "variante" CV, que nos permite elegir los mejores hiperparámetros y configuraciones de una forma sencilla.
Por ejemplo: la regresión Ridge tiene su variante RidgeCV que nos permite pasarle a la especificación del modelo una serie de hiperparámetros $\alpha$ a explorar (el hiperparámetro de shrinkage de la regularización $L2$). Sklearn entrenará cada modelo $K$ veces haciendo validación cruzada y nos devolverá el mejor modelo entrenado. Utilicémoslo con nuestro dataset de training de regresión (de años trabajados ~ salario anual):
from sklearn.linear_model import RidgeCV
ridge_cv = RidgeCV(alphas=[0.05, 0.1, 0.5, 1.0, 5.0, 10.0, 20.0], # Queremos probar con estos valores
fit_intercept=True,
cv=10 # 10-Fold Cross Validation
)
# Lo entrenamos. En el fondo, estamos entrenando
# 10 folds * 8 valores de alpha distintos = 80 modelos;
# así que tarda un poquito:
ridge_cv.fit(X=dataset_regresion_training[["años trabajados"]],
y=dataset_regresion_training["salario bruto"])
Modelos entrenados, y sklearn se ha "quedado" automáticamente con el mejor. El mejor de ellos ha sido el que tiene el hiperparámetro $\alpha$:
ridge_cv.alpha_
Ese es el hiperparámetro de los elegidos que minimiza el $MSE$ en validación cruzada.
Los modelos que tienen una variante CV consiguen hacer la validación cruzada de una forma muy eficiente, pero son poco flexibles (normalmente te dejan elegir los valores a explorar de uno o dos hiperparámetros, pero no de todo).
Si queremos trabajar más la validación cruzada, podemos utilizar directamente sklearn.cross_validation.KFold, que genera un iterador sobre los datos (para que usemos un bucle for, vayamos entrenando con distintos valores de los hiperparámetros y argumentos que queramos, y vayamos acumulando los resultados en una lista/diccionario por ejemplo).
Pero aún mejor: GridSearchCV nos permite hacer lo mismo de una forma más elegante y menos tediosa.
9.B.b. GridSearchCV¶
Grid Search es una técnica que nos permite construir una gran cantidad de modelos, entrenarlos y validarlos (con validación cruzada) de forma sencilla y expresiva.
El Grid Search está compuesto por:
- Un modelo (supervisado: clasificador o regresor)
- Un grid de hiperparámetros, donde definimos para cada hiperparámetro los valores que queremos probar
- Una técnica de validación cruzada
- Una métrica de rendimiento
Al igual que en los modelos con validación cruzada empotrada, Grid Search entrena los modelos con validación cruzada y probando, en este caso, las distintas combinaciones de hiperparámetros; y midiendo el rendimiento de cada modelo. El resultado final es el modelo que mejor rendimiento ha ofrecido.
GridSearchCV es la clase de sklearn que nos permite hacer Grid Search.
Vamos a probar en nuestro dataset de clasificación (calvos/no-calvos). Supongamos que queremos probar con K-Nearest Neighbors: a priori no sabemos qué valor de $K$ será mejor, así que podemos construir una Grid Search. El grid de hiperparámetros lo construimos con un diccionario, donde cada clave es un string con el nombre exacto de cada hiperparámetro (en el caso de sklearn.neighbors.KNeighborsClassifier es "n_neighbors"). Si el modelo toma más hiperparámetros/argumentos de los cuales no especificamos valores a explorar, los dejará constantes (y serán los que vienen por defecto en el modelo):
from sklearn.model_selection import GridSearchCV
knn_gs = GridSearchCV(estimator = KNeighborsClassifier(), # Instanciamos un KNN
param_grid = {"n_neighbors":[1,3,5,7,9,11]}, # Los valores de n_neighbors a probar, en lista
scoring="accuracy", # La métrica de rendimiento a utilizar
cv=5, # 5-Fold Cross Validation
verbose=3 # Esto hace que imprima "cosas informativas" durante la ejecución
)
knn_gs es ahora una instancia de GridSearchCV, que podemos tratar como un modelo normal y corriente (con sus métodos .fit() y .predict(), por ejemplo). Si aplicamos el método .fit():
knn_gs.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
pass
Nos ha impreso cómo ha ido entrenando y validando cada modelo con validación cruzada, el valor de n_neighbors que ha probado cada vez, el rendimiento obtenido y el tiempo que ha tardado en entrenar.
Podemos preguntarle a knn_gs cuál ha sido el mejor modelo en validación con el atributo .best_estimator_:
knn_gs.best_estimator_
Una vez entrenado, validado y elegido con el método .fit(), knn_gs utilizará el modelo ganador si aplicamos knn_gs.predict(), de forma que no es necesario que creemos otro KNeighborsClassifier con el número n_negihbors ganador para empezar a predecir con el mejor modelo. No obstante, lo recomiendo: utiliza GridSearchCV para conseguir la información de los hiperparámetros que mejor funcionan. Una vez que los conozcas, entrena el modelo de la forma convencional sin GridSearchCV; con los mejores hiperparámetros y argumentos.
Hagamos otro ejemplo de GridSearchCV, esta vez con árboles regresores sobre nuestro dataset de regresión (años trabajados ~ salario anual). Lo bueno es que los árboles toman muchos argumentos e hiperparámetros distintos; y GridSearchCV entrena todas las combinaciones posibles de todos ellos (que puede ser lento, pero es eficaz):
# Creamos el Grid Search de árboles:
arbol_gs = GridSearchCV(estimator=DecisionTreeRegressor(),
param_grid={"max_depth": [2,3,4,5,6],
"min_samples_split": [2,5,10,15],
"min_samples_leaf": [2,5,10,15]},
scoring="neg_mean_squared_error",
cv=5,
verbose=1 # Para que imprima menos información que antes
)
# Entrenamos, validamos y elegimos el mejor árbol:
arbol_gs.fit(X=dataset_regresion_training[["años trabajados"]],
y=dataset_regresion_training["salario bruto"])
pass
Vemos que a medida que metemos más hiperparámetros y más valores a explorar para cada uno de ellos, el número de modelos entrenados (y validados) aumenta exponencialmente. Veamos cuál es nuestro árbol ganador:
arbol_gs.best_estimator_
GridSearchCV nos ayda muchísimo a elegir los mejores hiperparámetros y argumentos para un modelo. Pero... ¿Y cuando queremos elegir entre varios modelos distintos (por ejemplo, entre el mejor árbol y el mejor SVM)?
Sencillo: resulta que cualquier GridSearchCV generado nos devuelve también el rendimiento del mejor modelo a través del atributo .best_score_:
arbol_gs.best_score_
Que en este caso es el $MSE$ (puesto que así lo definimos en arbol_gs al generarlo). Que no te despiste el signo menos que tiene delante el $MSE$: es por razones internas de sklearn. Básicamente, en la versión 0.18 (la que estamos utilizando aquí) se calcula el $-MSE$ (en negativo). Pero es exactamente lo mismo; solo que en vez de minimizar el $MSE$, se maximiza el $-MSE$, que es completamente equivalente.
El caso es que podemos:
- Construir varios
GridSearchCV: uno para cada modelo. Esto nos permitirá elegir los mejores hiperparámetros para cada modelo. - Meter todos los
GridSearchCVen una lista de Python normal y corriente. - Entrenar todos ellos con
.fit(), y obtener los.best_score_de cada modelo ganador. - Elegir el modelo con mejor
.best_score_
Vamos a ponerlo en práctica con nuestro problema de clasificación de la alopecia. Vamos a probar los siguientes modelos:
- Regresión logística
- K-NN
- SVM
- Random Forest
Buscando los mejores hiperparámetros para cada uno, y luego quedándonos con el mejor modelo. Manos a la obra:
# Definimos todos los modelos con su GridSearchCV:
reg_log_gs = GridSearchCV(estimator=LogisticRegression(),
param_grid={"penalty": ["l1","l2"], # Probamos con regularizaciones L1 y L2
"C": [0.5, 1.0, 1.5]}, # Y varios valores de la regularización C
scoring="roc_auc", # En vez de accuracy, vamos a usar el área bajo la curva ROC
cv=5)
knn_gs = GridSearchCV(estimator=KNeighborsClassifier(),
param_grid={"n_neighbors": [1,3,5,7,9,11]},
scoring="roc_auc", # Lo suyo es usar la misma métrica en todos los modelos que probemos
cv=5)
svm_gs = GridSearchCV(estimator=SVC(),
param_grid={"C": [0.5, 1.0, 1.5],
"kernel": ["linear","poly","rbf"],
"degree": [2,3,4]},
scoring="roc_auc",
cv=5)
random_f_gs = GridSearchCV(estimator=RandomForestClassifier(),
param_grid={"n_estimators": [5,10,20,50],
"max_depth": [2,3,4,5],
"min_samples_split": [3,4,5,10],
"min_samples_leaf": [3,4,5,10]},
scoring="roc_auc",
cv=5)
# Los metemos todos en una lista:
bag_de_modelos = [reg_log_gs, knn_gs, svm_gs, random_f_gs]
# Los entrenamos todos, validando y eligiendo el mejor de cada categoría
# (la mejor regresión logística, el mejor K-NN, el mejor SVM y el mejor Random Forest):
for modelo in bag_de_modelos:
modelo.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
y=dataset_clasificacion_training["label"])
# Creamos otra lista para meter los .best_score_:
lista_de_scores = []
# Y los metemos:
for modelo in bag_de_modelos:
lista_de_scores.append(modelo.best_score_)
Y los resultados son (recordamos que las listas de Python mantienen el orden de los elementos; así que el primer score es de la regresión logística, el segundo de K-NN, etc):
lista_de_scores
Ha ganado el tercer modelo, es decir, el SVM (SVC): es el que mayor área bajo la curva ROC ha dado. Vamos a ver con qué argumentos/hiperparámetros ha ganado:
bag_de_modelos[2].best_estimator_
Y de los generados, este es el modelo que querríamos utilizar para hacer nuestras predicciones a futuro.
La validación, así como la comparación y selección de modelos es la esencia del aprendizaje automático. Siempre buscamos un modelo capaz de conseguir buenos resultados con datos con los que no ha entrenado; y la validación cruzada junto con la búsqueda de hiperparámetros nos permite encontrar el modelo con el menor $E_{out}$ posible.
Por último, vamos a ver un elemento que cada vez está cobrando más importancia en la ciencia de datos: las pipelines.
10. Pipelines¶
Ya hemos visto las principales ramas de sklearn. Solo queda por comentar el tema de las pipelines.
Una pipeline no es más que una especie de diagrama que encadena distintos procesos. Normalmente, tenemos que seguir un determinado orden a la hora de utilizar aprendizaje automático: preprocesar y transformar los datos, seleccionar variables, entrenar modelos, validarlos y seleccionarlos. Las pipelines intentan simplificarnos la tarea: encadenan varias de estas fases, y generan un modelo portable, fácil de entrenar y entender.
Las pipelines son utilizadas en muchas bibliotecas de machine learning. Si alguna vez has usado herramientas visuales de aprendizaje automático y minería de datos (como SAS Enterprise Miner, SPSS Modeller, RapidMiner o KNIME), seguro habrás "arrastrado" distintos transformadores y modelos, los habrás configurado, y luego los habrás conectado entre sí. Pues una pipeline en sklearn es exactamente lo mismo.
Sklearn implementa sklearn.pipeline.Pipeline para construir pipelines.
Vamos a utilizar el dataset del precio de las casas de Boston. Querremos hacer lo siguiente, y por el siguiente orden:
- Seleccionar variables
- Puesto que las variables son de variopintos ámbitos y magnitudes, probablemente sea buena idea estandarizarlas (o igual no...)
- Entrenar, validar y seleccionar modelos
Vamos a hacerlo todo con una Pipeline:
from sklearn.pipeline import Pipeline
# Lo primero será generar un RFECV para seleccionar variables,
# que utilizará un árbol regresor sencillo para ver qué variables
# son las que debemos seleccionar:
arbol_seleccionar = DecisionTreeRegressor()
rfecv = RFECV(estimator=arbol_seleccionar)
# Luego utilizaremos un StandardScaler para estandarizar
# las variables:
scaler = StandardScaler()
# Luego probaremos diferentes random forest regresores:
random_forest_r = RandomForestRegressor()
# Ahora generamos la Pipeline, donde pasaremos una lista
# de las "steps" o pasos a seguir, a modo de tuplas:
# donde el primer elemento es un string con el "nombre"
# de cada cosa, y el segundo elemento es cada cosa en sí:
mi_pipeline = Pipeline(steps=[("seleccionador",rfecv), ("escalador",scaler), ("forest",random_forest_r)])
mi_pipeline
Hemos generado la pipeline. No obstante, vamos a querer hacer validación en todas las steps: en el modelo seleccionador de variables, en el escalador, y en el random forest. Para eso, podemos generar un diccionario que luego pasaremos a GridSearchCV, junto con la propia pipeline: utilizando los nombres personalizados que le hemos puesto a cada step, junto con la notación de poner doble underscore (__) y cada argumento/hiperparámetro, podemos definir los hiperparámetros que queremos probar de cada parte de la pipeline:
hiperparam_pipeline = {"seleccionador__cv":[5,10], # Si queremos explorar varios, los pasamos en lista
"seleccionador__scoring":["neg_mean_squared_error"], # Si no es así, lista con 1 elemento
"escalador__with_mean":[True, False], # Esta y la de abajo nos permiten explorar
"escalador__with_std":[True, False], # el usar el escalador, o que "no haga nada";
# que es lo mismo que no ponerlo; de esta forma,
# explorarmos si es mejor estandarizar o no
# en este caso. Si pones tanto with_mean como
# with_std en False en StandardScaler, es como
# si no estuviera, porque ni hace que la media
# de cada variable sea 0; ni que la desviación
# típica sea 1. También permite tocar la media
# y no la varianza, o al revés. Es decir: todas
# las combinaciones.
"forest__n_estimators":[10,20,40],
"forest__criterion":["mse"],
"forest__max_depth":[5,10,20],
"forest__min_samples_split":[2,3,5,10],
"forest__min_samples_leaf":[2,3,5]
}
Ahora podemos concebir mi_pipeline como "un gran modelo"; y por otro lado tenemos un diccionario con los hiperparámetros que queremos probar y validar en cada paso. Tal y como si fuera un modelo normal y corriente de sklearn, mi_pipeline tiene el método .fit(); que aplica .fit() en todos los elementos que tiene dentro.
Gracias a esto, podemos utilizar GridSearchCV para probar qué combinaciones son las mejores. Solo tenemos que pasar como modelo mi_pipeline, y como grid de hiperparámetros a explorar el diccionario hiperparametros_pipeline que hemos construido:
mi_pipeline_grid_search = GridSearchCV(estimator=mi_pipeline,
param_grid=hiperparam_pipeline,
scoring="neg_mean_squared_error",
cv=10,
n_jobs=4 # Para que use 4 CPUs a la vez; y entrene
# un modelo en cada una a la vez.
)
# Y le damos caña. Esto va a ser largo:
mi_pipeline_grid_search.fit(X=boston_training_set["features"],
y=boston_training_set["target","Precio"])
Vamos a ver cuál es el modelo (o mejor dicho: pipeline) ganador con .best_estimator_:
pipeline_ganadora = mi_pipeline_grid_search.best_estimator_
pipeline_ganadora
Vamos a ver los steps de la pipeline ganadora en más detalle con el atributo .steps, que devuelve una lista con los 3 procesos (en este caso):
pipeline_ganadora.steps
Y esa es la combinación ganadora. Ahora, lo ideal sería re-entrenar con los hiperparámetros que han resultado mejores, sobre todo el dataset de training (sin validación, porque ya la hemos hecho). Habiendo decidido que este es nuestro modelo (o pipeline) ganador, podríamos ver el rendimiento que da en el conjunto de testing; así tendríamos una idea de lo bien que debería predecir el modelo en el futuro.
Los modelos (así como los pre-procesamientos previos, selectores de variables...) deben ser re-entrenados con el tiempo. En el caso de que llegara una cantidad considerable de datos nuevos, lo ideal sería repetir el proceso desde 0. De esta forma "garantizamos" que nuestro trabajo sigue siendo útil a lo largo del tiempo.
11. Bibliografía y recursos adicionales¶
Dividiré esta sección en dos partes.
Por un lado, los libros que recomiendo relacionados con el aprendizaje automático y la minería de datos en general (no atados a Python, sino más a la matemática en la que se fundamenta cada modelo):
- The Elements of Statistical Learning: EL libro de aprendizaje automático. Trata una gran variedad de modelos y técnicas, explicadas en gran detalle. Obligatorio. Punto.
- Pattern Recognition and Machine Learning: también muy buen libro sobre aprendizaje automático. Eso sí: si eres muy frecuentista (y por tanto poco bayesiano) puedes acabar odiándolo. Recomendado a los amantes de la probabilidad y de Bayes, ya que intenta explicar todos los modelos a partir de ahí.
- Data Mining: un libro bastante reciente (2015). No solo entra en los modelos, sino que presta bastante atención a las cosas que los rodean (como el feature engineering y las métricas de rendimiento de modelos y de distancia).
- Bayesian Nonparametric Data Analysis: si quieres aprender más sobre cómo utilizar modelos bayesianos para el aprendizaje automático. Corto, pero intenso.
Y hay muchísimos más.
Luego ya dentro de Python hay buenos libros (como Building Machine Learning Systems with Python y Python Machine Learning) sobre machine learning con sklearn, además de otros más avanzados como Advanced Machine Learning with Python (sobre modelos menos convencionales y más avanzados) y Large Scale Machine Learning with Python (sobre machine learning a gran escala con Python) que son muy nuevos. Además, existen muchas bibliotecas más allá de sklearn útiles para el mundo del data science y machine learning:
- Bibliotecas de redes neuronales y Deep Learning: las comentamos anteriormente. Destaco Theano, TensorFlow y Keras.
- Bibliotecas para procesamiento de lenguaje natural: NLTK es una de las bibliotecas más famosas (si no la más) para hacer procesamiento de lenguaje natural. Ofrece muchas funciones y utilidades para el prepocesado de texto no estructurado (tokenizers, stemmers, computadores de n-grams, lemmatizers, part-of-speech taggers...) así como modelos de machine learning para utilizar una vez que el texto ha sido "traducido" a datos que los modelos pueden utilizar. No obstante, es frecuente utilizar NLTK para preprocesar el texto y utilizar luego sklearn para los modelos en sí (porque son más potentes). Genism es otra biblioteca muy buena; esencial para topic modelling y word embeddings (que incluye modelos como Word2vec). Los mejores libros para hacer procesamiento de lenguaje natural con Python son Natural Language Processing with Python (segunda edición en camino) y Python 3 Text Processing with NLTK 3 Cookbook.
- Bibliotecas para visión computacional y procesamiento de imágenes: OpenCV es la biblioteca más famosa para visión computacional: escrita en C++, incluye API para Python y hay buenos libros sobre el tema. Scikit-image también incluye una gran variedad de algoritmos para procesado de imágenes. A partir de ahí, es frecuente de nuevo utilizar sklearn o una biblioteca de deep learning para aprender de las imágenes.
- SymPy es una biblioteca de cálculo simbólico muy útil en ciertos escenarios relacionados con el aprendizaje automático, al igual que Theano (pero más sencilla).
- XGBoost es una biblioteca muy popular creada por la comunidad para formular y entrenar modelos de Gradient Boosted Trees. Es más escalable que la que viene en sklearn gracias a su muy ingeniosa implementación. Además, incluye una novedosa regularización que evita muy bien el sobreajuste en la práctica. Se han ganado muchos concursos de Kaggle utilizándola.
- Scipy incluye un módulo de optimización llamado
scipy.optimizeque puede ser muy útil cuando formulamos un modelo desde cero, y queremos tener un buen optimizador de la función. - StatsModels es una biblioteca estadística y econométrica que utiliza Numpy y Scipy. Personalmente no me gusta mucho, puesto que prefiero utilizar
scipy.statsdirectamente. Pero viene bien explorarla y ver de qué es capaz. - Spark (MLlib): Spark tiene una maravillosa API de Python que permite analizar y aplicar aprendizaje automático sobre datasets muy grandes de forma distribuida. La biblioteca MLlib de aprendizaje automático está inspirada en sklearn, hasta el punto que nos encontramos los métodos
.fit(),.predict()y.transform()por doquier. La Spark ML Pipeline es también muy similar a las Pipelines de sklearn. - Ibis pretende convertirse en un pandas distribuido, exponiendo una API muy similar a la de éste. Utiliza principalmente Impala como motor, y se prevé que en el futuro incluya modelos de machine learning. El ecosistema Blaze pretende hacer algo similar con Dask, pero menos atado al ecosistema Hadoop (de hecho, ya se pueden utilizar algunos modelos de sklearn en distribuido utilizando Dask).
Pero como digo siempre: ¡hay muchas más!
Profesor: Julio Antonio Soto